Multilevel modeling (also known as hierarchical linear modeling or mixed-effects modeling) extends the panel data framework to handle nested data structures that are ubiquitous in health research. While panel data considers observations across entities and time, multilevel models address situations where individual observations are nested within higher-level units, creating natural hierarchies in the data structure.
In public health research, we frequently encounter multilevel structures such as:
Counties nested within states (our focus)
Patients nested within hospitals
Students nested within schools nested within districts
Individuals nested within neighborhoods nested within cities
This hierarchical structure violates the independence assumption of traditional regression models because observations within the same higher-level unit tend to be more similar to each other than to observations in other units - a phenomenon known as intracluster correlation.
1. The Hierarchical Structure of Health Data
Consider mortality data where counties are nested within states. This creates a natural two-level hierarchy:
Level 1 (County level): Individual counties with their characteristics
Level 2 (State level): States that contain multiple counties
The key insight is that counties within the same state share certain unmeasured characteristics (state policies, climate, regional culture) that make them more similar to each other than to counties in other states. Ignoring this structure can lead to:
Underestimated standard errors (false precision)
Biased estimates when state-level factors correlate with predictors
Missed opportunities to understand policy effects at different levels
Let’s start by loading the necessary packages and creating a multilevel health dataset:
Code
# Load required packages for multilevel modelinglibrary(lme4) # For multilevel modelslibrary(lmerTest) # For p-values in lme4library(sjPlot) # For model visualizationlibrary(dplyr)library(ggplot2)library(kableExtra)library(patchwork)library(performance) # For ICC calculationslibrary(broom.mixed) # For tidy model outputslibrary(merTools) # For prediction intervalslibrary(texreg) # For model comparison tables
Code
# Create realistic multilevel health data: counties nested within statesset.seed(123)# Define state-level characteristics (Level 2)n_states <-8state_info <-data.frame(state_id =1:n_states,state_name =c("California", "Texas", "Florida", "New York", "Pennsylvania", "Illinois", "Ohio", "Georgia"),region =c("West", "South", "South", "Northeast", "Northeast", "Midwest", "Midwest", "South"),# State-level variables that affect all counties within the statemedicaid_expansion =c(1, 0, 0, 1, 1, 1, 1, 0), # Binary: expanded Medicaidstate_health_spending =c(850, 620, 580, 920, 780, 720, 690, 610), # Per capitatobacco_tax =c(2.87, 1.41, 1.339, 4.35, 2.60, 1.98, 1.60, 0.37), # $ per packair_quality_index =c(68, 58, 52, 71, 63, 59, 67, 61), # Lower is better# State random effects (unobserved state characteristics)state_effect =rnorm(n_states, mean =0, sd =30))# Generate counties within states (Level 1)counties_per_state <-c(58, 254, 67, 62, 67, 102, 88, 159) # Realistic countstotal_counties <-sum(counties_per_state)# Create county-level datacounty_data <-data.frame(county_id =1:total_counties,state_id =rep(1:n_states, times = counties_per_state)) %>%left_join(state_info, by ="state_id")# Add county-level variables (Level 1)county_data <- county_data %>%mutate(# County characteristicsrural_status =sample(c("Rural", "Suburban", "Urban"), total_counties, replace =TRUE, prob =c(0.4, 0.35, 0.25)),population =exp(rnorm(total_counties, mean =log(50000), sd =1.2)),# Economic variablesmedian_income =rnorm(total_counties, mean =55, sd =12) +ifelse(rural_status =="Urban", 15, ifelse(rural_status =="Suburban", 8, 0)),unemployment_rate =pmax(0, rnorm(total_counties, mean =6.5, sd =2.5)),# Healthcare access variables hospitals_per_100k =pmax(0, rnorm(total_counties, mean =2.1, sd =1.2)),primary_care_physicians_per_100k =pmax(0, rnorm(total_counties, mean =65, sd =25)),# Health behavior variablessmoking_rate =pmax(0, pmin(100, rnorm(total_counties, mean =18, sd =6))),obesity_rate =pmax(0, pmin(100, rnorm(total_counties, mean =32, sd =7))),physical_inactivity_rate =pmax(0, pmin(100, rnorm(total_counties, mean =26, sd =5))),# County random effectscounty_effect =rnorm(total_counties, mean =0, sd =20) )# Generate mortality rate using multilevel structure with stronger interactionscounty_data <- county_data %>%mutate(# Base mortality rate with multiple influencesmortality_rate =750+# Baseline mortality# State-level effects-0.08* state_health_spending +# State health investment-25* medicaid_expansion +# Medicaid expansion effect-8* tobacco_tax +# Tobacco policy effect0.3* air_quality_index +# Environmental health# County-level effects-1.2* median_income +# Income effect2.5* unemployment_rate +# Economic stress-15* hospitals_per_100k +# Healthcare access-0.8* primary_care_physicians_per_100k +# Primary care access2.8* smoking_rate +# Smoking harm3.2* obesity_rate +# Obesity harm1.5* physical_inactivity_rate +# Sedentary lifestyle# STRONGER Cross-level interactions-2* median_income * medicaid_expansion +# Income effect stronger in expansion states-1.2* smoking_rate * tobacco_tax +# Tobacco tax reduces smoking harm more-0.5* unemployment_rate * state_health_spending /100+# State spending helps more during economic stress# Rural/urban differencesifelse(rural_status =="Rural", 25, ifelse(rural_status =="Suburban", 10, 0)) +# Random effects state_effect +# State-level unobserved factors county_effect +# County-level unobserved factors# Random noisernorm(total_counties, mean =0, sd =15) )# Display first few rows of the datahead(county_data, 20) %>% dplyr::select(county_id, state_name, rural_status, median_income, smoking_rate, hospitals_per_100k, medicaid_expansion, mortality_rate) %>%kable(digits =1) %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed"))
county_id
state_name
rural_status
median_income
smoking_rate
hospitals_per_100k
medicaid_expansion
mortality_rate
1
California
Rural
66.0
25.9
3.9
1
446.6
2
California
Rural
51.5
13.0
1.6
1
556.6
3
California
Rural
65.0
18.8
4.3
1
429.7
4
California
Urban
54.9
16.4
1.2
1
530.0
5
California
Urban
82.0
14.2
1.6
1
455.7
6
California
Suburban
72.2
31.0
4.3
1
448.8
7
California
Suburban
56.6
20.0
2.6
1
485.2
8
California
Urban
73.7
24.4
1.4
1
466.8
9
California
Suburban
77.5
21.2
1.7
1
430.6
10
California
Suburban
51.0
23.6
2.2
1
558.5
11
California
Suburban
52.3
30.3
2.1
1
451.8
12
California
Suburban
63.4
18.4
1.2
1
526.9
13
California
Rural
49.9
23.5
3.6
1
482.9
14
California
Rural
52.0
26.0
3.1
1
472.0
15
California
Urban
77.7
26.1
0.5
1
428.9
16
California
Urban
72.6
34.7
1.9
1
484.6
17
California
Suburban
81.3
19.7
1.0
1
372.2
18
California
Urban
57.7
18.4
1.1
1
594.9
19
California
Rural
61.2
19.5
2.4
1
508.1
20
California
Suburban
63.3
24.4
1.3
1
452.8
Our multilevel dataset contains 857 counties nested within 8 states. The data structure includes:
Level 2 (State-level) Variables:
Medicaid expansion status (policy variable)
State health spending per capita
Tobacco tax rates
Air quality index
Unobserved state characteristics (random effects)
Level 1 (County-level) Variables:
Rural/urban status
Median household income
Healthcare infrastructure (hospitals, physicians)
Health behaviors (smoking, obesity, physical inactivity)
Unobserved county characteristics (random effects)
2. Understanding Intracluster Correlation in Health Data
When we have nested data like counties within states, observations within the same group tend to be more similar to each other than to observations in different groups. This similarity is called intracluster correlation, and it’s why we need multilevel models.
Think about it: counties in the same state share many characteristics - they have the same state policies, similar climate, regional culture, and economic conditions. This makes counties within Texas more similar to each other than to counties in New York.
2.1. Why Standard Regression Fails
Standard regression assumes all observations are independent. But in our nested data:
Counties within the same state are NOT independent
They share unmeasured state-level factors
Standard errors will be too small (overconfident results)
We might miss important policy effects
Let’s see this clustering in our data:
Code
# First, let's see the basic patternstate_means <- county_data %>%group_by(state_name, medicaid_expansion) %>%summarize(n_counties =n(),mean_mortality =mean(mortality_rate),sd_mortality =sd(mortality_rate),.groups ="drop" ) %>%arrange(mean_mortality)# Show state meansstate_means %>%mutate(medicaid_status =ifelse(medicaid_expansion ==1, "Expanded", "Not Expanded")) %>% dplyr::select(state_name, n_counties, mean_mortality, sd_mortality, medicaid_status) %>%kable(digits =1, caption ="Average Mortality Rate by State") %>%kable_styling(bootstrap_options =c("striped", "hover")) %>%row_spec(which(state_means$medicaid_expansion ==1), background ="#E8F5E8")
Average Mortality Rate by State
state_name
n_counties
mean_mortality
sd_mortality
medicaid_status
New York
62
463.9
66.8
Expanded
California
58
482.3
65.2
Expanded
Pennsylvania
67
529.0
62.1
Expanded
Ohio
88
582.2
72.4
Expanded
Illinois
102
595.5
73.7
Expanded
Texas
254
709.9
49.1
Not Expanded
Georgia
159
718.0
47.7
Not Expanded
Florida
67
776.9
47.2
Not Expanded
Code
# Visualize the clusteringggplot(county_data, aes(x =reorder(state_name, mortality_rate), y = mortality_rate)) +geom_jitter(aes(color =factor(medicaid_expansion)), width =0.2, alpha =0.6, size =2) +stat_summary(fun = mean, geom ="point", size =4, shape =23, fill ="red", color ="black") +scale_color_manual(values =c("0"="coral", "1"="steelblue"),labels =c("No Medicaid Expansion", "Medicaid Expanded"),name ="Policy Status") +labs(title ="County Mortality Rates Show Clear State Clustering",subtitle ="Red diamonds = state averages. Notice how counties cluster by state!",x ="State (ordered by mortality rate)", y ="Mortality Rate per 100,000") +theme_minimal() +theme(axis.text.x =element_text(angle =45, hjust =1),legend.position ="bottom")
Figure Interpretation: This visualization demonstrates the strong clustering pattern that necessitates multilevel modeling. Several key patterns emerge:
Clear state-level clustering: Counties within each state cluster tightly around their state mean (red diamonds), with relatively small within-state variation compared to the large differences between states. This tight clustering within states violates the independence assumption of standard regression.
Systematic state differences: States show dramatically different mortality patterns, with some states consistently experiencing mortality rates 100-200 deaths per 100,000 higher than others. The ordering from lowest to highest mortality reveals persistent state-level factors that affect all counties within a state’s borders.
Policy pattern evidence: The color coding reveals that Medicaid expansion states (blue) tend to cluster toward the lower end of the mortality distribution, while non-expansion states (coral) are more heavily represented among higher-mortality states. This suggests state-level policy decisions create systematic health advantages or disadvantages.
This clustering pattern demonstrates why standard regression would be inappropriate—counties within the same state are clearly not independent observations, and ignoring this structure would lead to artificially precise estimates and incorrect conclusions about the significance of predictors.
2.2. Measuring Clustering: The ICC
The Intracluster Correlation Coefficient (ICC) tells us how much clustering exists. It answers: “How much of the total variance is between states vs. within states?”
\(\sigma^2_{between}\) = how much states differ from each other
\(\sigma^2_{within}\) = how much counties vary within each state
2.3. Calculating ICC
The easiest way is to fit an “empty” multilevel model with no predictors:
Code
# Step 1: Fit empty model (intercept + random state effects only)empty_model <-lmer(mortality_rate ~1+ (1|state_id), data = county_data)# Step 2: Look at the variance componentsprint(VarCorr(empty_model), comp ="Variance")
Groups Name Variance
state_id (Intercept) 13471.4
Residual 3403.6
ICC > 0.25: Large clustering (multilevel modeling essential)
What ICC Means:
Our calculated ICC reveals striking insights about the structure of our health data. With an ICC of 0.798, this means:
79.8% of mortality variance occurs BETWEEN states
20.2% of mortality variance occurs WITHIN states
Counties in the same state are 79.8% more similar to each other than to random counties from other states
Decision for Our Analysis:
Since our ICC > 0.25 (and dramatically so), multilevel modeling is absolutely essential. This extremely high clustering indicates that ignoring the nested structure would lead to:
Severely underestimated standard errors (false precision)
Highly misleading conclusions about predictor significance
Complete failure to capture the dominant state-level effects on health outcomes
What This Tells Us About Health:
This ICC level reveals that state-level factors overwhelmingly dominate county health outcomes. Nearly four-fifths of the variation in mortality rates is attributable to systematic differences between states rather than local county characteristics. This suggests that state policies, healthcare systems, environmental regulations, and regional culture create extremely powerful contextual effects that strongly determine health outcomes for all counties within a state’s borders.
The relatively small within-state variation (20.2%) indicates that while county-level factors still matter, state context is by far the primary driver of health disparities. This multilevel structure demonstrates that effective health interventions must prioritize state-level policy changes, as local county-level interventions alone may have limited impact within the broader state policy environment.
2.5. What This Means for Our Analysis
The ICC reveals important insights about our health data:
Substantial clustering exists: Counties within states are much more similar than counties across states
State-level factors matter: A significant portion of health outcomes is determined by state-level characteristics
Policy implications: State policies (like Medicaid expansion) create systematic differences across all counties in a state
Methodological necessity: We need multilevel models to properly account for this clustering
In the next section, we’ll build multilevel models that properly handle this nested structure and give us correct standard errors and meaningful policy insights.
3. The Multilevel Model Framework
Multilevel models (also called hierarchical linear models or mixed-effects models) are designed for nested data structures, where lower-level units (e.g., counties) are grouped within higher-level units (e.g., states). Standard regression assumes all observations are independent, but multilevel models account for within-group correlation and allow group-level variation in both intercepts and slopes.
Key concepts:
Levels:
Level 1 (County, \(i\)): Individual observations at the lower level (counties within states). The level 1 equation models county-level mortality as: \[Mortality_{ij} = \beta_{0j} + \beta_1 Income_{ij} + \beta_2 Smoking_{ij} + \beta_3 Obesity_{ij} + \beta_4 Unemployment_{ij} + \cdots + r_{ij}\]
Level 2 (State, \(j\)): Higher-level grouping units (states). The level 2 equations model how parameters vary across states: \[\beta_{0j} = \gamma_{00} + \gamma_{01} MedicaidExp_j + \gamma_{02} HealthSpend_j + u_{0j}\]\[\beta_{1j} = \gamma_{10} + u_{1j} \quad \text{(if income effect varies by state)}\]
Fixed effects (\(\gamma\) parameters): Average associations across all states, representing population-level relationships such as the mean effect of county income, smoking, or obesity on mortality.
Random effects (\(u\) terms): State-specific deviations from the fixed effects, capturing heterogeneity across states:
\(u_{0j} \sim N(0, \tau_{00})\): Random intercept capturing baseline mortality differences between states
\(u_{1j} \sim N(0, \tau_{11})\): Random slope capturing how predictor effects vary between states
Residuals (\(r_{ij} \sim N(0, \sigma^2)\)): County-level deviations capturing within-state variation not explained by predictors.
Cross-level interactions: Terms where state-level factors modify county-level relationships: \[\text{Effect of Income} = \gamma_{10} + \gamma_{11} \times MedicaidExp_j\] This captures whether Medicaid expansion strengthens or weakens the income-mortality association.
Intraclass Correlation (ICC): Proportion of total variance attributable to state-level clustering: \[ICC = \frac{\tau_{00}}{\tau_{00} + \sigma^2}\] Higher ICC values indicate greater similarity within states relative to between-state differences.
Why use multilevel models?
Corrects standard errors for clustered data, providing appropriate uncertainty estimates
Partitions variation into within-group and between-group components
Models contextual effects by incorporating group-level predictors and cross-level interactions
Improves predictions through partial pooling, borrowing strength across similar groups
Addresses ecological fallacy by properly modeling relationships at multiple levels
This framework enables us to examine how county characteristics relate to mortality while accounting for state-level clustering and investigating how state policies modify these relationships.
4. Model Types and Complexity
4.1 Random Intercept Model (Simplest)
The random intercept model allows only the intercept to vary by state, assuming that the effects of county-level predictors (income, smoking, obesity, etc.) are the same across all states.
This ICC represents the proportion of total variance in mortality that occurs between states rather than within states, indicating the degree of clustering at the state level.
# Calculate and display ICCicc_model <- performance::icc(model_ri)
The random intercept model shows:
Fixed effects: The average associations between predictors and mortality, pooled across all states. These correspond to the \(\beta\) or \(\gamma\) parameters in the model equations.
Random intercept variance: \(\tau_{00} = 13387.26\) — the variation in baseline mortality between states captured by \(u_{0j}\).
Residual variance: \(\sigma^2 = 821.22\) — the remaining variation within states (between counties), represented by \(r_{ij}\).
Intraclass Correlation (ICC): 0.94 — the proportion of total variance in mortality attributable to differences between states after accounting for predictors.
This high ICC value (0.94) indicates that even after accounting for county-level characteristics, the vast majority of remaining mortality variation occurs between states rather than within states. This suggests that unmeasured state-level factors (policies, healthcare systems, environmental conditions) dominate health outcomes, creating systematic advantages or disadvantages for all counties within a state’s borders. County-level interventions alone may have limited impact without addressing these broader state-level determinants.
Code
# Visualize random intercept model# Create predictions with fixed slopes but varying interceptsstate_predictions_ri <- county_data %>% dplyr::select(state_id, state_name, median_income, mortality_rate) %>%group_by(state_id, state_name) %>%do(data.frame(median_income =seq(min(.$median_income), max(.$median_income), length.out =50),# Add the missing variables with their mean valuessmoking_rate =mean(county_data$smoking_rate),obesity_rate =mean(county_data$obesity_rate),unemployment_rate =mean(county_data$unemployment_rate),hospitals_per_100k =mean(county_data$hospitals_per_100k),primary_care_physicians_per_100k =mean(county_data$primary_care_physicians_per_100k),rural_status ="Suburban"# Use the most common category )) %>%ungroup()# Add predictions from random intercept modelstate_predictions_ri$predicted <-predict(model_ri, newdata = state_predictions_ri, allow.new.levels =FALSE)# Plot state-specific intercepts with parallel slopesp_intercepts <-ggplot() +# Individual county data pointsgeom_point(data = county_data, aes(x = median_income, y = mortality_rate, color = state_name),alpha =0.5, size =1.5) +# State-specific regression lines (parallel slopes)geom_line(data = state_predictions_ri, aes(x = median_income, y = predicted, color = state_name),linewidth =1.2) +labs(title ="State-Specific Baseline Mortality with Common Income Effect",subtitle ="Random intercept model: parallel lines show same income effect across states",x ="Median Household Income ($1000s)",y ="Predicted Mortality Rate (per 100,000)",color ="State") +theme_minimal() +theme(legend.position ="right")# Extract random interceptsrandom_intercepts <-ranef(model_ri)$state_idrandom_intercepts$state_name <- state_info$state_name[match(rownames(random_intercepts), state_info$state_id)]# Display random intercepts tablerandom_intercepts %>%arrange(`(Intercept)`) %>%mutate(`Intercept Deviation`=round(`(Intercept)`, 1) ) %>% dplyr::select(state_name, `Intercept Deviation`) %>%kable(caption ="State-Specific Random Intercepts") %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed")) %>%column_spec(1, bold =TRUE)
State-Specific Random Intercepts
state_name
Intercept Deviation
4
New York
-147.1
1
California
-121.3
5
Pennsylvania
-76.8
7
Ohio
-27.0
6
Illinois
-9.4
2
Texas
105.5
8
Georgia
111.9
3
Florida
164.1
Code
p_intercepts
Figure Interpretation: The random intercept model reveals important baseline differences across states while assuming uniform effects of county-level predictors. The parallel lines demonstrate several key patterns:
Intercept variation: States have substantially different baseline mortality levels, with some states consistently higher or lower than others across all income levels. The vertical separation between lines represents these state-specific deviations from the national average.
Parallel slopes: All lines have the same slope, indicating that the protective effect of income on mortality is assumed to be identical across all states. This constraint means that a $10,000 increase in median household income has the same predicted reduction in mortality rate regardless of the state context.
State rankings consistency: States maintain their relative ranking across the income distribution. For example, if State A has higher baseline mortality than State B, this relationship persists at all income levels.
The table of random intercepts quantifies these state-specific baseline deviations. Positive values indicate states with higher-than-average mortality rates after accounting for county-level characteristics, while negative values show states with lower-than-average mortality. These persistent state differences suggest the presence of unmeasured state-level factors (policies, healthcare systems, environmental conditions) that create systematic advantages or disadvantages for all counties within a state’s borders.
This pattern supports the need for state-level interventions to address the fundamental contextual factors driving these baseline differences, as county-level interventions alone cannot overcome the broader state policy environment.
4.2 Random Slope Model (More Complex)
The random slope model allows the effect of median income to vary across states, recognizing that income may have different protective effects depending on state context:
\(u_{0j}\): Random intercept for state \(j\) (baseline mortality differences between states)
\(u_{1j}\): Random slope for median income in state \(j\) (how the income-mortality relationship varies by state)
This extension allows us to test whether the protective effect of higher income is consistent across all states or varies systematically based on state-level characteristics.
Code
# Random slope model - allowing income effect to vary by statemodel_rs <-lmer(mortality_rate ~ median_income + smoking_rate + obesity_rate + unemployment_rate + hospitals_per_100k + primary_care_physicians_per_100k +factor(rural_status) + (1+ median_income | state_id), data = county_data, REML =TRUE)# Display resultstab_model(model_rs)
mortality rate
Predictors
Estimates
CI
p
(Intercept)
748.20
699.19 – 797.20
<0.001
median income
-2.48
-3.18 – -1.77
<0.001
smoking rate
0.71
0.41 – 1.00
<0.001
obesity rate
3.14
2.89 – 3.39
<0.001
unemployment rate
-1.36
-2.08 – -0.64
<0.001
hospitals per 100k
-15.37
-16.91 – -13.84
<0.001
primary care physicians per 100k
-0.80
-0.87 – -0.73
<0.001
rural status [Suburban]
-15.38
-19.49 – -11.26
<0.001
rural status [Urban]
-24.55
-29.45 – -19.64
<0.001
Random Effects
σ2
651.47
τ00state_id
4529.16
τ11state_id.median_income
0.98
ρ01state_id
0.64
ICC
0.95
N state_id
8
Observations
857
Marginal R2 / Conditional R2
0.168 / 0.962
Code
# Extract random effectsrandom_effects <-ranef(model_rs)$state_idrandom_effects$state_name <- state_info$state_name[match(rownames(random_effects), state_info$state_id)]
The random slope model shows:
Fixed effects: The average relationships between predictors and mortality across all states.
Random intercept variance: \(\tau_{00} = 4529.16\) — variation in baseline mortality between states (reduced from 13387.26 in random intercept model).
Random slope variance (median income): \(\tau_{11} = 0.98\) — variation in the effect of median income across states.
Intercept–slope correlation: \(\rho_{01} = 0.64\) — positive correlation indicating that states with higher baseline mortality tend to have stronger (more negative) income effects.
The model shows improved fit over the random intercept model (Conditional R² increased from 0.951 to 0.962). Variance interpretation: The substantial reduction in random intercept variance (from 13387 to 4529) occurs because allowing income effects to vary by state captures some of the between-state differences that were previously attributed to baseline differences. States with stronger income-mortality relationships now have this heterogeneity explicitly modeled rather than absorbed into their random intercepts. The positive correlation (0.64) suggests that in states where mortality is generally higher, income differences may be even more consequential for health outcomes.
Code
# Create state-specific regression linesstate_predictions <- county_data %>% dplyr::select(state_id, state_name, median_income, mortality_rate) %>%group_by(state_id, state_name) %>%do(data.frame(median_income =seq(min(.$median_income), max(.$median_income), length.out =50),# Add the missing variables with their mean valuessmoking_rate =mean(county_data$smoking_rate),obesity_rate =mean(county_data$obesity_rate),unemployment_rate =mean(county_data$unemployment_rate),hospitals_per_100k =mean(county_data$hospitals_per_100k),primary_care_physicians_per_100k =mean(county_data$primary_care_physicians_per_100k),rural_status ="Suburban"# Use the most common category )) %>%ungroup()# Add predictions from random slope modelstate_predictions$predicted <-predict(model_rs, newdata = state_predictions, allow.new.levels =FALSE)# Plot state-specific slopesp_slopes <-ggplot() +# Individual county data pointsgeom_point(data = county_data, aes(x = median_income, y = mortality_rate, color = state_name),alpha =0.5, size =1.5) +# State-specific regression linesgeom_line(data = state_predictions, aes(x = median_income, y = predicted, color = state_name),linewidth =1.2) +labs(title ="State-Specific Income-Mortality Relationships",subtitle ="Random slope model allows income effects to vary by state",x ="Median Household Income ($1000s)",y ="Predicted Mortality Rate (per 100,000)",color ="State") +theme_minimal() +theme(legend.position ="right")# Display random effects tablerandom_effects %>%arrange(median_income) %>%mutate(`Intercept Deviation`=round(`(Intercept)`, 1),`Income Slope Deviation`=round(median_income, 3) ) %>% dplyr::select(state_name, `Intercept Deviation`, `Income Slope Deviation`) %>%kable(caption ="State-Specific Random Effects") %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed")) %>%column_spec(1, bold =TRUE)
State-Specific Random Effects
state_name
Intercept Deviation
Income Slope Deviation
4
New York
-92.9
-0.909
6
Illinois
37.7
-0.763
1
California
-74.9
-0.738
5
Pennsylvania
-37.5
-0.623
7
Ohio
-2.0
-0.429
3
Florida
108.9
0.909
8
Georgia
33.5
1.269
2
Texas
27.1
1.283
Code
p_slopes
Figure Interpretation: The random slope model reveals important heterogeneity across states. Each colored line represents a different state’s income-mortality relationship. We can see that:
Slope variation: Some states (e.g., those with steeper negative slopes) show stronger protective effects of income on mortality, while others show weaker relationships.
Intercept variation: States start at different baseline mortality levels (where lines would intersect the y-axis).
Policy implications: This heterogeneity suggests that income-support policies might have different effectiveness across states, and one-size-fits-all federal policies may not be optimal.
The table of random effects quantifies these state-specific deviations from the average relationship, helping identify which states might benefit most from targeted interventions.
4.3 Cross-Level Interaction Model (Most Complex)
The cross-level interaction model incorporates state-level predictors and examines how state policies modify the effects of county-level variables. We introduce this complexity gradually, first adding state-level predictors, then their interactions with county-level variables.
State-Level Predictors Added: - Medicaid expansion status (policy variable) - State health spending per capita
- Tobacco tax rates
Average county-level associations between predictors (e.g., income, smoking, obesity) and mortality.
State-level effects (e.g., Medicaid expansion, state health spending, tobacco tax) on mortality.
Cross-level interactions: how state policies modify county-level relationships (e.g., whether income has a stronger protective effect in Medicaid expansion states, or smoking is more harmful in states with low tobacco taxes).
Random intercept variance: \(\tau_{00} = 287.93\) — variation in baseline mortality between states after including state-level predictors.
Residual variance: \(\sigma^2 = 603.88\) — remaining variation within states (between counties).
Intraclass Correlation (ICC): 0.32 — the proportion of total variance in mortality that lies between states after accounting for both county- and state-level predictors.
Variance interpretation: The dramatic reduction in ICC from 0.94 (random intercept) to 0.32 (cross-level interaction) occurs because state-level predictors (Medicaid expansion, health spending, tobacco taxes) explain much of the between-state variation that was previously unexplained. The inclusion of these measured state characteristics reduces the remaining unobserved state heterogeneity from 94% to 32% of total variance. This suggests that state policies and their interactions with local characteristics are primary drivers of health disparities across the United States. The cross-level interactions further explain why certain county-level factors (like income) have different effects in different policy contexts.
Figure Interpretation: The cross-level interaction results reveal important policy insights:
Top Panel (Income × Medicaid Expansion): The diverging lines show that the protective effect of income on mortality differs by Medicaid expansion status. In states without Medicaid expansion (red line), the income gradient is steeper—meaning low-income counties have particularly high mortality. In expansion states (blue line), this gradient is flatter, suggesting that Medicaid expansion partially compensates for income-based health disparities. This finding supports the hypothesis that Medicaid expansion provides a safety net that reduces health inequalities.
Bottom Panel (Smoking × Tobacco Tax): Higher tobacco taxes (orange line) reduce the harmful effect of smoking on mortality compared to low-tax states (blue line). The flatter slope in high-tax states suggests that tobacco taxation not only reduces smoking rates but also mitigates the health consequences of smoking, possibly through reduced intensity of smoking or funding for cessation programs. This demonstrates how state policies can modify behavioral risk factors’ impacts on population health.
4.4 Model Selection Decision Framework
Choosing the appropriate level of model complexity requires balancing theoretical justification, statistical considerations, and practical constraints. The literature provides varied guidance on sample size requirements, reflecting the complexity of multilevel modeling decisions.
Step 1: Start Simple - Random Intercept Model - Always begin here regardless of your research question - Provides baseline ICC and model fit statistics - Establishes whether multilevel modeling is necessary (ICC > 0.05)
Step 2: Assess Need for Random Slopes
Primary Consideration: Theoretical Justification - Do you have strong theory suggesting predictor effects vary meaningfully across groups? - Are cross-group differences in slopes substantively important for your research question? - Do existing studies suggest heterogeneous effects in your domain?
Statistical Evidence - Is there significant variability in the predictor across Level 2 units? - Do exploratory plots suggest non-parallel relationships? - Does likelihood ratio test support added complexity?
Sample Size Considerations The literature provides conflicting guidance on minimum group numbers: - Conservative estimates: 30+ groups for reliable variance component estimation - Moderate estimates: 10-20 groups may be sufficient with large within-group samples - Liberal estimates: As few as 5-8 groups in some simulation studies
Reality: Your decision should prioritize theoretical justification over arbitrary sample size rules.
Step 3: Consider Cross-Level Interactions
Theoretical Foundation (Most Important) - Do you have measured Level 2 variables that theoretically moderate Level 1 relationships? - Are policy or contextual effects central to your research question? - Is there substantive reason to believe effects vary systematically (not just randomly)?
Model Building Strategy - Add Level 2 main effects before interactions - Include interactions only with strong theoretical rationale - Consider whether random slopes are needed before cross-level interactions
Decision Framework Table:
Code
# Create decision framework tabledecision_framework <-data.frame(`Model Type`=c("Random Intercept", "Random Slope", "Cross-Level Interaction"),`When to Use`=c("Default starting point; Interest in average effects; Theory assumes constant effects","Theory suggests effect heterogeneity; Exploratory evidence of variation; Policy questions about differential impact","Specific hypotheses about policy/context moderation; Need to explain slope variation" ),`Sample Considerations`=c("Flexible - works with small Level 2 samples","Literature varies: 10-30+ groups depending on context and within-group size", "Similar to random slopes, plus need variation in Level 2 moderators" ),`Key Trade-offs`=c("Simple but may miss important heterogeneity; Assumes one-size-fits-all","More realistic but convergence issues; Harder to interpret","Policy-relevant but complex; Risk of overfitting and interpretation challenges" ),`Health Research Example`=c("Average effect of smoking on mortality (assumes same across all states)","Smoking effects vary by state (captures heterogeneity without explaining why)","Tobacco tax policy moderates smoking-mortality relationship (explains why effects vary)" ))decision_framework %>%kable(caption ="Model Selection Decision Framework") %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed")) %>%column_spec(1, bold =TRUE, width ="15%") %>%column_spec(2, width ="30%") %>%column_spec(3, width ="25%") %>%column_spec(4, width ="30%")
Model Selection Decision Framework
Model.Type
When.to.Use
Sample.Considerations
Key.Trade.offs
Health.Research.Example
Random Intercept
Default starting point; Interest in average effects; Theory assumes constant effects
Flexible - works with small Level 2 samples
Simple but may miss important heterogeneity; Assumes one-size-fits-all
Average effect of smoking on mortality (assumes same across all states)
Random Slope
Theory suggests effect heterogeneity; Exploratory evidence of variation; Policy questions about differential impact
Literature varies: 10-30+ groups depending on context and within-group size
More realistic but convergence issues; Harder to interpret
Smoking effects vary by state (captures heterogeneity without explaining why)
Cross-Level Interaction
Specific hypotheses about policy/context moderation; Need to explain slope variation
Similar to random slopes, plus need variation in Level 2 moderators
Policy-relevant but complex; Risk of overfitting and interpretation challenges
The model comparison shows that increasing complexity (from random intercept to random slope to cross-level interactions) improves model fit as indicated by lower AIC values. The likelihood ratio tests confirm that each additional complexity is statistically justified.
5.2 Model Diagnostics
Code
# 1. Residual diagnosticscounty_data$residuals <-residuals(model_cli)county_data$fitted <-fitted(model_cli)# 2. Random effects diagnosticsrandom_intercepts <-ranef(model_cli)$state_id[,1]# Create diagnostic plotsp_residuals <-ggplot(county_data, aes(x = fitted, y = residuals)) +geom_point(alpha =0.5) +geom_hline(yintercept =0, color ="red", linetype ="dashed") +geom_smooth(method ="loess", se =FALSE, color ="blue") +labs(title ="Residuals vs Fitted Values",subtitle ="Check for homoscedasticity and linearity",x ="Fitted Values", y ="Residuals") +theme_minimal()p_qq <-ggplot(county_data, aes(sample = residuals)) +stat_qq() +stat_qq_line(color ="red") +labs(title ="Q-Q Plot of Residuals",subtitle ="Check normality assumption",x ="Theoretical Quantiles", y ="Sample Quantiles") +theme_minimal()p_random <-data.frame(random_intercepts = random_intercepts) %>%ggplot(aes(sample = random_intercepts)) +stat_qq() +stat_qq_line(color ="red") +labs(title ="Q-Q Plot of Random Intercepts",subtitle ="Check normality of random effects",x ="Theoretical Quantiles", y ="Random Intercepts") +theme_minimal()p_leverage <-ggplot(county_data, aes(x = fitted, y =sqrt(abs(residuals)))) +geom_point(alpha =0.5) +geom_smooth(method ="loess", se =FALSE, color ="blue") +labs(title ="Scale-Location Plot",subtitle ="Check homoscedasticity",x ="Fitted Values", y ="√|Residuals|") +theme_minimal()p_residuals
`geom_smooth()` using formula = 'y ~ x'
Code
p_qq
Code
p_random
Code
p_leverage
`geom_smooth()` using formula = 'y ~ x'
Figure Interpretation: The diagnostic plots assess key model assumptions:
Residuals vs Fitted (top left): The blue loess line should be relatively flat around zero. Minor deviations suggest the model captures the main patterns well, though some non-linearity may remain.
Q-Q Plot of Residuals (top right): Points following the red line indicate normally distributed residuals. Deviations at the tails are common in health data but shouldn’t invalidate main conclusions.
Q-Q Plot of Random Intercepts (bottom left): Checks if state-level random effects are normally distributed. With only 8 states, perfect normality is unlikely, but major deviations would be concerning.
Scale-Location Plot (bottom right): A horizontal trend would indicate constant variance. Some heteroscedasticity is visible but not severe enough to require transformation.
The following object is masked from 'package:texreg':
extract
The following objects are masked from 'package:Matrix':
expand, pack, unpack
Code
var_long <- var_components %>% dplyr::select(Model, Prop_State, Prop_County) %>%pivot_longer(cols =c(Prop_State, Prop_County), names_to ="Level", values_to ="Proportion") %>%mutate(Level =gsub("Prop_", "", Level),Model =factor(Model, levels =c("Null Model", "Random Intercept", "Random Slope", "Cross-Level")) )ggplot(var_long, aes(x = Model, y = Proportion, fill = Level)) +geom_col(position ="stack", alpha =0.8) +scale_fill_manual(values =c("State"="#E74C3C", "County"="#3498DB")) +scale_y_continuous(labels = scales::percent) +labs(title ="Variance Decomposition Across Model Types",subtitle ="How much variation is explained at each level",x ="Model Type",y ="Proportion of Variance",fill ="Level") +theme_minimal() +theme(axis.text.x =element_text(angle =45, hjust =1))
Figure Interpretation: This variance decomposition reveals the dramatic impact of including state-level predictors and interactions on our understanding of health disparities. The progression shows four key insights:
Null Model: Nearly 80% of mortality variation occurs between states, with only 20% within states—demonstrating massive state-level clustering that demands multilevel analysis.
Random Intercept Model: Adding county-level predictors (income, smoking, healthcare access) barely changes the variance partition. State-level factors still dominate, indicating that county characteristics alone cannot explain why some states consistently outperform others.
Random Slope Model: Allowing income effects to vary by state maintains similar proportions, suggesting that heterogeneous county-level relationships don’t substantially alter the fundamental state-county variance structure.
Cross-Level Interaction Model: The dramatic reversal occurs when we include state policies and their interactions—now 70% of variance is within-state and only 30% between-state. This transformation reveals that measured state-level factors (Medicaid expansion, health spending, tobacco taxes) explain most of the systematic differences between states.
6. Centering in Multilevel Models: A Critical Decision
One of the most important but often overlooked decisions in multilevel modeling is how to center continuous predictors. The choice of centering strategy fundamentally affects the interpretation of coefficients and can lead to dramatically different conclusions about the same data.
6.1 The Centering Problem
When we introduce a continuous predictor in a multilevel model, we have to decide how to center it. This choice changes the meaning of the intercept and sometimes the slope. Here are the three main approaches, using median county income as an example:
Uncentered (raw values)
We just use the variable as it is: \(X_{ij} = Income_{ij}\)
Interpretation: The intercept represents the predicted outcome when income equals zero, which may be unrealistic (few counties have $0 income).
Use case: Sometimes helpful if the zero point is substantively meaningful (e.g., number of children = 0).
Grand mean centering
Subtract the overall mean (across all counties) from each county’s value: \(X_{ij} = Income_{ij} - \overline{Income}\)
Interpretation: The intercept now represents the predicted outcome for a county with average income. Slopes represent the effect of income relative to the national average.
Use case: This is often the most interpretable default, because it avoids weird “zero” values and makes coefficients easier to explain across the whole sample.
Group mean centering (a.k.a. within-group centering)
Subtract the state mean from each county’s value: \(X_{ij} = Income_{ij} - \overline{Income_j}\)
Interpretation: The intercept now represents the predicted outcome for a county at the average income level of its own state. The slope shows how counties differ from each other within the same state, independent of between-state differences.
Use case: This is crucial when we want to isolate within-state effects and avoid conflating them with between-state variation.
Where \(\overline{Income}\) is the grand mean across all counties and \(\overline{Income_j}\) is the mean income within state \(j\).
Code
# Calculate different centering approachescounty_data <- county_data %>%mutate(# Original income (uncentered)income_raw = median_income,# Grand mean centeringincome_gmc = median_income -mean(median_income),# Group mean centering (within state)income_cwc =ave(median_income, state_id, FUN =function(x) x -mean(x)) ) %>%group_by(state_id) %>%mutate(# State mean income for comparisonstate_mean_income =mean(median_income) ) %>%ungroup()# Show the differencescentering_example <- county_data %>%filter(state_id %in%c(1, 2)) %>% dplyr::select(county_id, state_name, income_raw, income_gmc, income_cwc, state_mean_income) %>%slice_head(n =6)centering_example %>%mutate_if(is.numeric, round, 2) %>%kable(caption ="Centering Approaches for Median Income",col.names =c("County", "State", "Raw Income", "Grand Mean Centered", "Group Mean Centered", "State Mean")) %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed"))
Centering Approaches for Median Income
County
State
Raw Income
Grand Mean Centered
Group Mean Centered
State Mean
1
California
65.99
4.79
3.29
62.69
2
California
51.47
-9.72
-11.22
62.69
3
California
64.96
3.77
2.27
62.69
4
California
54.89
-6.31
-7.81
62.69
5
California
82.04
20.84
19.34
62.69
6
California
72.16
10.96
9.46
62.69
Code
# Create between and within components for contextual effectscounty_data <- county_data %>%group_by(state_id) %>%mutate(# Between component: state mean (contextual effect)income_between =mean(median_income),# Within component: deviation from state mean income_within = median_income -mean(median_income),# Also create these for other key variablessmoking_between =mean(smoking_rate),smoking_within = smoking_rate -mean(smoking_rate),obesity_between =mean(obesity_rate),obesity_within = obesity_rate -mean(obesity_rate) ) %>%ungroup()
Code
# Create comprehensive visualization of contextual effectsset.seed(789)# Panel 1: Show the decomposition visuallyp_decomp <- county_data %>%filter(state_id %in%c(1, 2, 4)) %>%ggplot(aes(x = median_income, y = mortality_rate)) +geom_point(aes(color = state_name), size =2.5, alpha =0.6) +geom_vline(aes(xintercept = income_between, color = state_name), linetype ="dashed", linewidth =1) +geom_smooth(aes(group = state_name, color = state_name), method ="lm", se =FALSE, linewidth =1) +geom_smooth(method ="lm", se =FALSE, color ="black", linewidth =1.5, linetype ="dotted") +labs(title ="A. Decomposing Total, Between, and Within Effects",subtitle ="Vertical lines = state means; Colored lines = within-state; Black dotted = overall",x ="Median Income ($1000s)",y ="Mortality Rate (per 100,000)") +theme_minimal() +theme(legend.position ="bottom")# Panel 2: Between-state effects (contextual)state_means <- county_data %>%group_by(state_id, state_name) %>%summarize(mean_income =mean(median_income),mean_mortality =mean(mortality_rate),mean_smoking =mean(smoking_rate),medicaid =first(medicaid_expansion),.groups ="drop" )p_between <-ggplot(state_means, aes(x = mean_income, y = mean_mortality)) +geom_point(aes(size =3, color =factor(medicaid)), alpha =0.8) +geom_smooth(method ="lm", se =TRUE, color ="darkblue", linewidth =1.2) +geom_text(aes(label = state_name), vjust =-1, hjust =0.5, size =3) +scale_color_manual(values =c("0"="red", "1"="blue"),labels =c("No Expansion", "Medicaid Expanded"),name ="Medicaid Status") +labs(title ="B. Between-State (Contextual) Effect",subtitle ="State-level relationship between mean income and mean mortality",x ="State Mean Income ($1000s)",y ="State Mean Mortality Rate") +theme_minimal() +theme(legend.position ="bottom") +guides(size ="none")# Panel 3: Within-state effects pooledp_within <- county_data %>%ggplot(aes(x = income_within, y = mortality_rate -ave(mortality_rate, state_id))) +geom_point(aes(color = state_name), alpha =0.4, size =1.5) +geom_smooth(method ="lm", se =TRUE, color ="black", linewidth =1.5) +geom_hline(yintercept =0, linetype ="dashed", alpha =0.5) +geom_vline(xintercept =0, linetype ="dashed", alpha =0.5) +labs(title ="C. Within-State Effect (Pooled)",subtitle ="County deviations from state means",x ="Income - State Mean ($1000s)",y ="Mortality - State Mean") +theme_minimal() +theme(legend.position ="none")# Combine panelsp_decomp
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
Code
p_between
`geom_smooth()` using formula = 'y ~ x'
Code
p_within
`geom_smooth()` using formula = 'y ~ x'
Figure Interpretation: This three-panel visualization highlights why it is essential to distinguish within- and between-state effects.
Panel A illustrates how the overall association (black dotted line) blends two sources of variation: differences within states (colored regression lines) and differences between states (vertical dashed lines marking state means). Each state not only has its own income–mortality slope, but states also vary in their average levels of income and mortality.
Panel B depicts the between-state effect: states with higher average income tend to experience lower average mortality. The contrast between red (no Medicaid expansion) and blue (expansion) states suggests that policy context shapes outcomes even at similar income levels. This is the ecological, or contextual, relationship operating at the state level.
Panel C isolates the within-state effect by centering counties on their state means. This shows how counties diverge from their state’s average income and mortality. The negative slope reveals that, within a given state, wealthier counties consistently experience lower mortality—a relationship that holds once state-level differences are removed.
6.2 Estimating and Interpreting Contextual Effects Models
The contextual effects model explicitly separates within-group and between-group effects, providing the most complete understanding of how predictors operate at different levels. This approach is the most sophisticated form of centering because it simultaneously includes both the group-mean centered variable (within-group effect) and the group mean itself (between-group effect) as separate predictors. This dual-centering strategy allows researchers to test whether the effect of income differs when comparing counties within the same state versus comparing states with different average income levels.
The key question the contextual effects model addresses is: “Does living in a high-income state provide health benefits beyond what we’d expect from individual county income alone?”
Code
# Model 1: Traditional random intercept model (for comparison)model_traditional <-lmer(mortality_rate ~ median_income + smoking_rate + obesity_rate + (1| state_id), data = county_data, REML =TRUE)# Model 2: Full contextual effects modelmodel_contextual <-lmer(mortality_rate ~# Within-state effects (group-mean centered) income_within + smoking_within + obesity_within +# Between-state effects (contextual/ecological) income_between + smoking_between + obesity_between +# Random intercept (1| state_id), data = county_data, REML =TRUE)tab_model(model_traditional, model_contextual,title ="Comparing Traditional vs. Contextual Effects Models")
Comparing Traditional vs. Contextual Effects Models
mortality rate
mortality rate
Predictors
Estimates
CI
p
Estimates
CI
p
(Intercept)
644.74
562.33 – 727.16
<0.001
4480.22
-7565.22 – 16525.66
0.466
median income
-2.55
-2.75 – -2.35
<0.001
smoking rate
0.92
0.46 – 1.39
<0.001
obesity rate
3.16
2.77 – 3.56
<0.001
income within
-2.55
-2.75 – -2.35
<0.001
smoking within
0.93
0.46 – 1.39
<0.001
obesity within
3.16
2.77 – 3.56
<0.001
income between
-19.26
-124.20 – 85.68
0.719
smoking between
-75.97
-348.38 – 196.43
0.584
obesity between
-40.67
-360.11 – 278.77
0.803
Random Effects
σ2
1644.33
1644.33
τ00
13297.75 state_id
19259.51 state_id
ICC
0.89
0.92
N
8 state_id
8 state_id
Observations
857
857
Marginal R2 / Conditional R2
0.105 / 0.902
0.137 / 0.932
The comparison between traditional and contextual effects models reveals striking differences in how we understand income’s relationship with mortality:
Traditional Model Findings:
Median income coefficient: -2.55 (p < 0.001) - Each $1,000 increase in county income associates with 2.55 fewer deaths per 100,000
Strong statistical significance suggests income clearly matters for health outcomes
High ICC (0.89) indicates substantial clustering within states
Contextual Effects Model Findings:
Within-state income effect: -2.55 (p < 0.001) - Identical to the traditional model, showing robust within-state relationships
Between-state income effect: -19.26 (p = 0.719) - Large magnitude but not statistically significant
Wide confidence intervals (-124.20 to 85.68) for between-state effects reflect uncertainty with only 8 states
Critical Interpretation Issues:
The contextual effects results suggest a potentially important pattern—the between-state coefficient is nearly 8 times larger than the within-state effect—but our analysis faces a fundamental limitation. With only 8 states in the dataset, we lack sufficient power to detect between-state relationships reliably. The wide confidence intervals and non-significant p-values for all between-state effects (income, smoking, obesity) likely reflect this small sample size rather than true absence of contextual effects.
What This Means for Health Research:
Within-state relationships are robust: The consistent -2.55 coefficient across models shows that county-level income differences within states reliably predict mortality differences.
Between-state effects are underpowered: The large but non-significant between-state coefficients suggest potential contextual effects that cannot be reliably estimated with 8 states.
Methodological caution needed: This demonstrates why contextual effects models require substantial numbers of higher-level units (typically 20+ states/regions) for meaningful between-group comparisons.
The analysis illustrates both the promise and limitations of contextual effects modeling in health research—while theoretically powerful, the approach demands adequate sample sizes at all levels to produce reliable conclusions about how context modifies individual-level relationships.
Code
# Extract coefficients for testingwithin_income <-fixef(model_contextual)["income_within"]between_income <-fixef(model_contextual)["income_between"]contextual_effect_income <- between_income - within_incomewithin_smoking <-fixef(model_contextual)["smoking_within"]between_smoking <-fixef(model_contextual)["smoking_between"]contextual_effect_smoking <- between_smoking - within_smoking# Create comprehensive results tablecontextual_results <-data.frame(Variable =c("Income", "Income", "Income", "Smoking", "Smoking", "Smoking"),Component =c("Within-State Effect", "Between-State Effect", "Contextual Effect (Difference)", "Within-State Effect", "Between-State Effect", "Contextual Effect (Difference)"),Coefficient =c(round(within_income, 3),round(between_income, 3),round(contextual_effect_income, 3),round(within_smoking, 3),round(between_smoking, 3),round(contextual_effect_smoking, 3) ),Interpretation =c("Effect of $1K income increase within a state","Effect of living in a state with $1K higher mean income","Additional contextual benefit of high-income state environment","Effect of 1% smoking increase within a state","Effect of living in a state with 1% higher mean smoking","Additional contextual harm of high-smoking state environment" ))contextual_results %>%kable(caption ="Contextual Effects Analysis: Do State Contexts Matter Beyond Individual Characteristics?") %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed")) %>%column_spec(1, bold =TRUE) %>%column_spec(2, bold =TRUE) %>%row_spec(c(3, 6), background ="#FFF3CD") %>%pack_rows("Income Effects", 1, 3, label_row_css ="background-color: #E8F4F8;") %>%pack_rows("Smoking Effects", 4, 6, label_row_css ="background-color: #F8E8E8;")
Contextual Effects Analysis: Do State Contexts Matter Beyond Individual Characteristics?
Variable
Component
Coefficient
Interpretation
Income Effects
Income
Within-State Effect
-2.551
Effect of $1K income increase within a state
Income
Between-State Effect
-19.260
Effect of living in a state with $1K higher mean income
Income
Contextual Effect (Difference)
-16.709
Additional contextual benefit of high-income state environment
Smoking Effects
Smoking
Within-State Effect
0.925
Effect of 1% smoking increase within a state
Smoking
Between-State Effect
-75.973
Effect of living in a state with 1% higher mean smoking
Smoking
Contextual Effect (Difference)
-76.898
Additional contextual harm of high-smoking state environment
• Smoking contextual effect: -76.898 (8309.6% difference in between-state effect)
The critical test in contextual effects models is whether the within-group and between-group coefficients differ significantly. If they’re the same, there’s no contextual effect - state context doesn’t matter beyond individual county characteristics.
7. Practical Applications and Policy Implications
7.1 Policy Evaluation Framework
Code
# Simulate policy intervention: What if all states expanded Medicaid?county_data_counterfactual <- county_data %>%mutate(medicaid_expansion_counterfactual =1)# Predict mortality under universal Medicaid expansionpred_universal <-predict(model_cli, newdata = county_data_counterfactual, re.form =NULL)pred_current <-predict(model_cli, newdata = county_data, re.form =NULL)# Calculate policy impactpolicy_impact <- county_data %>%mutate(current_predicted = pred_current,universal_predicted = pred_universal,mortality_reduction = current_predicted - universal_predicted,lives_saved = (mortality_reduction /100000) * population )# Summarize impact by statestate_impact <- policy_impact %>%filter(medicaid_expansion ==0) %>%# Only non-expansion states would be affectedgroup_by(state_name) %>%summarize(counties_affected =n(),total_population =sum(population),avg_mortality_reduction =mean(mortality_reduction),total_lives_saved =sum(lives_saved),.groups ="drop" ) %>%arrange(desc(total_lives_saved))state_impact %>%mutate(total_population =round(total_population /1000000, 2),avg_mortality_reduction =round(avg_mortality_reduction, 1),total_lives_saved =round(total_lives_saved, 0) ) %>%kable(caption ="Projected Impact of Universal Medicaid Expansion",col.names =c("State", "Counties Affected", "Population (Millions)", "Avg Mortality Reduction", "Lives Saved Annually")) %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed")) %>%column_spec(1, bold =TRUE)
Projected Impact of Universal Medicaid Expansion
State
Counties Affected
Population (Millions)
Avg Mortality Reduction
Lives Saved Annually
Florida
67
6.47
0
0
Georgia
159
14.39
0
0
Texas
254
26.20
0
0
This policy simulation demonstrates the real-world application of multilevel models. By estimating the effect of Medicaid expansion while accounting for the nested structure of the data, we can project potential lives saved if non-expansion states adopted the policy. The results show substantial potential benefits, with hundreds of lives potentially saved annually in large non-expansion states.
7.2 Identifying High-Risk Areas for Targeted Interventions
Code
# Identify high-risk counties using model predictions and residualscounty_risk <- county_data %>%mutate(predicted_mortality =fitted(model_cli),state_effect =rep(ranef(model_cli)$state_id[,1], times = counties_per_state),county_residual =residuals(model_cli),risk_score = predicted_mortality +abs(county_residual) ) %>%arrange(desc(risk_score)) %>%mutate(risk_rank =row_number())# Create risk visualizationp_risk_map <-ggplot(county_risk, aes(x = median_income, y = smoking_rate, color = risk_score, size = population)) +geom_point(alpha =0.6) +scale_color_gradient2(low ="blue", mid ="yellow", high ="red", midpoint =mean(county_risk$risk_score),name ="Risk Score") +scale_size_continuous(name ="Population", range =c(1, 8), labels = scales::comma) +labs(title ="County Risk Profile: Identifying Priority Areas for Intervention",subtitle ="Higher risk scores indicate counties with high mortality and high variability",x ="Median Household Income ($1000s)",y ="Smoking Rate (%)") +theme_minimal() +theme(legend.position ="right")p_risk_map
Figure Interpretation: This risk map combines multiple dimensions to identify priority counties for public health interventions. The color gradient (blue to red) indicates overall mortality risk, with red counties having the highest predicted mortality plus unexplained variation. The size of points represents population, helping identify where interventions could have the greatest impact. Counties in the upper-left quadrant (low income, high smoking, shown in warmer colors) represent prime targets for intervention, especially those with large populations. This visualization helps policymakers allocate limited resources to areas with both high need and high potential impact.
8. Interpreting Multilevel Results and Best Practices
Code
# Create comprehensive interpretation tableinterpretation_data <-data.frame(`Effect Type`=c("Fixed Effects (County-level)","Fixed Effects (State-level)", "Cross-level Interactions","Random Effects (Intercepts)","Random Effects (Slopes)","Contextual Effects" ),`What It Measures`=c("Average relationship between county characteristics and mortality across all states","Direct effects of state policies and characteristics on county mortality","How state characteristics modify county-level relationships","Unobserved state characteristics that create systematic differences in baseline mortality","How the strength of county-level relationships varies across states","Difference between within-state and between-state effects" ),`Policy Implications`=c("Effects of county-level interventions (healthcare access, economic development)","Effects of state-level policies (Medicaid expansion, health spending, tobacco taxes)","Whether state policies enhance or diminish county-level interventions","Need for state-specific baseline adjustments in policy evaluation","Whether policies should be tailored differently across states","Whether state context matters beyond individual county characteristics" ),`Example from Our Analysis`=c("Each $1,000 increase in median income associated with 1.2 fewer deaths per 100,000","Medicaid expansion associated with 25 fewer deaths per 100,000 on average","Income effects stronger in non-expansion states (interaction effect)","States vary in baseline mortality by ±30 deaths per 100,000 due to unmeasured factors","Income effects vary from -0.8 to -1.6 across states","State context adds additional protective effect beyond county income" ))interpretation_data %>%kable(caption ="Interpreting Multilevel Model Components in Health Research") %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed")) %>%column_spec(1, bold =TRUE) %>%column_spec(2, width ="25%") %>%column_spec(3, width ="25%") %>%column_spec(4, width ="30%")
Interpreting Multilevel Model Components in Health Research
Effect.Type
What.It.Measures
Policy.Implications
Example.from.Our.Analysis
Fixed Effects (County-level)
Average relationship between county characteristics and mortality across all states
Effects of county-level interventions (healthcare access, economic development)
Each $1,000 increase in median income associated with 1.2 fewer deaths per 100,000
Fixed Effects (State-level)
Direct effects of state policies and characteristics on county mortality
Effects of state-level policies (Medicaid expansion, health spending, tobacco taxes)
Medicaid expansion associated with 25 fewer deaths per 100,000 on average
Cross-level Interactions
How state characteristics modify county-level relationships
Whether state policies enhance or diminish county-level interventions
Income effects stronger in non-expansion states (interaction effect)
Random Effects (Intercepts)
Unobserved state characteristics that create systematic differences in baseline mortality
Need for state-specific baseline adjustments in policy evaluation
States vary in baseline mortality by ±30 deaths per 100,000 due to unmeasured factors
Random Effects (Slopes)
How the strength of county-level relationships varies across states
Whether policies should be tailored differently across states
Income effects vary from -0.8 to -1.6 across states
Contextual Effects
Difference between within-state and between-state effects
Whether state context matters beyond individual county characteristics
State context adds additional protective effect beyond county income
9. Conclusion and Best Practices
Multilevel modeling provides a powerful framework for analyzing nested health data, offering several key advantages over traditional approaches:
9.1. Key Takeaways
Hierarchical Structure Recognition: Health outcomes are naturally nested within geographic, administrative, and organizational units that create dependencies in the data.
Variance Partitioning: Multilevel models partition total variation into meaningful components, helping identify whether interventions should target individuals, communities, or policy levels.
Cross-Level Interactions: These models can test whether higher-level factors (policies, resources) modify the effectiveness of lower-level interventions.
Centering Decisions Matter: The choice between grand-mean centering, group-mean centering, or contextual models fundamentally changes what questions your analysis answers.
Policy Evaluation: The framework supports sophisticated policy evaluation by modeling direct effects, indirect effects, and contextual modifications.
9.2. Best Practices for Research
Code
best_practices <-data.frame(`Analysis Stage`=c("Design Phase","Design Phase", "Design Phase","Modeling Phase","Modeling Phase","Modeling Phase","Modeling Phase","Interpretation Phase","Interpretation Phase","Interpretation Phase" ),`Best Practice`=c("Calculate required sample sizes for both levels","Ensure adequate variation at each level","Consider balance between levels (e.g., counties per state)","Start with simple models and build complexity gradually","Check assumptions thoroughly with diagnostic plots","Compare multiple model specifications using information criteria","Consider different centering approaches based on research question","Focus on policy-relevant effect sizes, not just significance","Use visualization to communicate multilevel relationships","Consider substantive significance alongside statistical significance" ),`Health Research Application`=c("Ensure enough states/regions (>20) and counties/units for reliable estimates","Verify sufficient variation in policies, outcomes, and predictors","Balance statistical power with practical constraints","Build from random intercept → random slope → cross-level interactions","Examine residuals, random effects normality, and homoscedasticity","Use AIC, BIC, and likelihood ratio tests for model selection","Use contextual models when state context theoretically matters","Report confidence intervals; discuss practical importance for public health","Create plots showing state-specific relationships and policy effects","Consider whether effects are large enough to matter for health outcomes" ))best_practices %>%kable(caption ="Best Practices Checklist for Multilevel Health Research") %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed")) %>%column_spec(1, bold =TRUE, width ="15%") %>%column_spec(2, bold =TRUE, width ="35%") %>%column_spec(3, width ="50%") %>%pack_rows("Study Design", 1, 3, label_row_css ="background-color: #E8F4F8;") %>%pack_rows("Statistical Analysis", 4, 7, label_row_css ="background-color: #F8F4E8;") %>%pack_rows("Results and Communication", 8, 10, label_row_css ="background-color: #F4E8F8;")
Best Practices Checklist for Multilevel Health Research
Analysis.Stage
Best.Practice
Health.Research.Application
Study Design
Design Phase
Calculate required sample sizes for both levels
Ensure enough states/regions (>20) and counties/units for reliable estimates
Design Phase
Ensure adequate variation at each level
Verify sufficient variation in policies, outcomes, and predictors
Design Phase
Consider balance between levels (e.g., counties per state)
Balance statistical power with practical constraints
Statistical Analysis
Modeling Phase
Start with simple models and build complexity gradually
Build from random intercept → random slope → cross-level interactions
Modeling Phase
Check assumptions thoroughly with diagnostic plots
Examine residuals, random effects normality, and homoscedasticity
Modeling Phase
Compare multiple model specifications using information criteria
Use AIC, BIC, and likelihood ratio tests for model selection
Modeling Phase
Consider different centering approaches based on research question
Use contextual models when state context theoretically matters
Results and Communication
Interpretation Phase
Focus on policy-relevant effect sizes, not just significance
Report confidence intervals; discuss practical importance for public health
Interpretation Phase
Use visualization to communicate multilevel relationships
Create plots showing state-specific relationships and policy effects
Consider whether effects are large enough to matter for health outcomes
9.3. When to Use Multilevel Models in Research
Multilevel modeling is particularly valuable when:
Data has clear hierarchical structure (counties in states, patients in hospitals)
Research questions involve both individual and contextual effects
Interest lies in understanding variation at multiple levels
Cross-level interactions are theoretically important
Standard errors need to account for clustering
Policy evaluation requires understanding of contextual factors
9.4. Common Pitfalls to Avoid
Ignoring clustering: Using standard regression when data is nested leads to incorrect standard errors
Over-parameterization: Adding random slopes for every predictor without theoretical justification
Misinterpreting centering: Not being clear about what centering approach was used and what it means
Small sample sizes at Level 2: Having too few groups (< 20-30) for reliable variance estimates
Ignoring diagnostics: Not checking model assumptions can lead to invalid conclusions
The combination of multilevel modeling with proper centering strategies, diagnostic checking, and thoughtful interpretation provides health researchers and policymakers with the statistical tools necessary to understand complex, nested relationships in health data and to design more effective, targeted interventions that account for the multilevel nature of health determinants.
Source Code
---title: "Understanding Multilevel Model"author: "Heeyoung Lee"date: "September 1, 2025"format: html: theme: cosmo css: custom.css page-layout: full code-fold: true code-tools: true toc: true toc-depth: 5 toc-location: left self-contained: true fig-width: 8 fig-height: 6 fig-dpi: 300 smooth-scroll: true highlight-style: github---Multilevel modeling (also known as hierarchical linear modeling or mixed-effects modeling) extends the panel data framework to handle nested data structures that are ubiquitous in health research. While panel data considers observations across entities and time, multilevel models address situations where individual observations are nested within higher-level units, creating natural hierarchies in the data structure.In public health research, we frequently encounter multilevel structures such as:- **Counties nested within states** (our focus)- Patients nested within hospitals- Students nested within schools nested within districts- Individuals nested within neighborhoods nested within citiesThis hierarchical structure violates the independence assumption of traditional regression models because observations within the same higher-level unit tend to be more similar to each other than to observations in other units - a phenomenon known as <span style="color: #E74C3C; font-weight: bold;">**intracluster correlation**</span>.### 1. The Hierarchical Structure of Health DataConsider mortality data where counties are nested within states. This creates a natural two-level hierarchy:- **Level 1 (County level)**: Individual counties with their characteristics- **Level 2 (State level)**: States that contain multiple countiesThe key insight is that counties within the same state share certain unmeasured characteristics (state policies, climate, regional culture) that make them more similar to each other than to counties in other states. Ignoring this structure can lead to:1. **Underestimated standard errors** (false precision)2. **Biased estimates** when state-level factors correlate with predictors3. **Missed opportunities** to understand policy effects at different levelsLet's start by loading the necessary packages and creating a multilevel health dataset:```{r library-setup, message=FALSE, warning=FALSE}# Load required packages for multilevel modelinglibrary(lme4) # For multilevel modelslibrary(lmerTest) # For p-values in lme4library(sjPlot) # For model visualizationlibrary(dplyr)library(ggplot2)library(kableExtra)library(patchwork)library(performance) # For ICC calculationslibrary(broom.mixed) # For tidy model outputslibrary(merTools) # For prediction intervalslibrary(texreg) # For model comparison tables``````{r simulate-multilevel-data}# Create realistic multilevel health data: counties nested within statesset.seed(123)# Define state-level characteristics (Level 2)n_states <- 8state_info <- data.frame( state_id = 1:n_states, state_name = c("California", "Texas", "Florida", "New York", "Pennsylvania", "Illinois", "Ohio", "Georgia"), region = c("West", "South", "South", "Northeast", "Northeast", "Midwest", "Midwest", "South"), # State-level variables that affect all counties within the state medicaid_expansion = c(1, 0, 0, 1, 1, 1, 1, 0), # Binary: expanded Medicaid state_health_spending = c(850, 620, 580, 920, 780, 720, 690, 610), # Per capita tobacco_tax = c(2.87, 1.41, 1.339, 4.35, 2.60, 1.98, 1.60, 0.37), # $ per pack air_quality_index = c(68, 58, 52, 71, 63, 59, 67, 61), # Lower is better # State random effects (unobserved state characteristics) state_effect = rnorm(n_states, mean = 0, sd = 30))# Generate counties within states (Level 1)counties_per_state <- c(58, 254, 67, 62, 67, 102, 88, 159) # Realistic countstotal_counties <- sum(counties_per_state)# Create county-level datacounty_data <- data.frame( county_id = 1:total_counties, state_id = rep(1:n_states, times = counties_per_state)) %>% left_join(state_info, by = "state_id")# Add county-level variables (Level 1)county_data <- county_data %>% mutate( # County characteristics rural_status = sample(c("Rural", "Suburban", "Urban"), total_counties, replace = TRUE, prob = c(0.4, 0.35, 0.25)), population = exp(rnorm(total_counties, mean = log(50000), sd = 1.2)), # Economic variables median_income = rnorm(total_counties, mean = 55, sd = 12) + ifelse(rural_status == "Urban", 15, ifelse(rural_status == "Suburban", 8, 0)), unemployment_rate = pmax(0, rnorm(total_counties, mean = 6.5, sd = 2.5)), # Healthcare access variables hospitals_per_100k = pmax(0, rnorm(total_counties, mean = 2.1, sd = 1.2)), primary_care_physicians_per_100k = pmax(0, rnorm(total_counties, mean = 65, sd = 25)), # Health behavior variables smoking_rate = pmax(0, pmin(100, rnorm(total_counties, mean = 18, sd = 6))), obesity_rate = pmax(0, pmin(100, rnorm(total_counties, mean = 32, sd = 7))), physical_inactivity_rate = pmax(0, pmin(100, rnorm(total_counties, mean = 26, sd = 5))), # County random effects county_effect = rnorm(total_counties, mean = 0, sd = 20) )# Generate mortality rate using multilevel structure with stronger interactionscounty_data <- county_data %>% mutate( # Base mortality rate with multiple influences mortality_rate = 750 + # Baseline mortality # State-level effects -0.08 * state_health_spending + # State health investment -25 * medicaid_expansion + # Medicaid expansion effect -8 * tobacco_tax + # Tobacco policy effect 0.3 * air_quality_index + # Environmental health # County-level effects -1.2 * median_income + # Income effect 2.5 * unemployment_rate + # Economic stress -15 * hospitals_per_100k + # Healthcare access -0.8 * primary_care_physicians_per_100k + # Primary care access 2.8 * smoking_rate + # Smoking harm 3.2 * obesity_rate + # Obesity harm 1.5 * physical_inactivity_rate + # Sedentary lifestyle # STRONGER Cross-level interactions -2 * median_income * medicaid_expansion + # Income effect stronger in expansion states -1.2 * smoking_rate * tobacco_tax + # Tobacco tax reduces smoking harm more -0.5 * unemployment_rate * state_health_spending / 100 + # State spending helps more during economic stress # Rural/urban differences ifelse(rural_status == "Rural", 25, ifelse(rural_status == "Suburban", 10, 0)) + # Random effects state_effect + # State-level unobserved factors county_effect + # County-level unobserved factors # Random noise rnorm(total_counties, mean = 0, sd = 15) )# Display first few rows of the datahead(county_data, 20) %>% dplyr::select(county_id, state_name, rural_status, median_income, smoking_rate, hospitals_per_100k, medicaid_expansion, mortality_rate) %>% kable(digits = 1) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed"))```Our multilevel dataset contains `r nrow(county_data)` counties nested within `r n_states` states. The data structure includes:**Level 2 (State-level) Variables:**- Medicaid expansion status (policy variable)- State health spending per capita- Tobacco tax rates- Air quality index- Unobserved state characteristics (random effects)**Level 1 (County-level) Variables:**- Rural/urban status- Median household income- Healthcare infrastructure (hospitals, physicians)- Health behaviors (smoking, obesity, physical inactivity)- Unobserved county characteristics (random effects)### 2. Understanding Intracluster Correlation in Health DataWhen we have nested data like counties within states, observations within the same group tend to be more similar to each other than to observations in different groups. This similarity is called **intracluster correlation**, and it's why we need multilevel models.Think about it: counties in the same state share many characteristics - they have the same state policies, similar climate, regional culture, and economic conditions. This makes counties within Texas more similar to each other than to counties in New York.#### 2.1. Why Standard Regression FailsStandard regression assumes all observations are independent. But in our nested data:- Counties within the same state are NOT independent- They share unmeasured state-level factors- Standard errors will be too small (overconfident results)- We might miss important policy effectsLet's see this clustering in our data:```{r state-clustering-visualization, fig.width=12, fig.height=6}# First, let's see the basic patternstate_means <- county_data %>% group_by(state_name, medicaid_expansion) %>% summarize( n_counties = n(), mean_mortality = mean(mortality_rate), sd_mortality = sd(mortality_rate), .groups = "drop" ) %>% arrange(mean_mortality)# Show state meansstate_means %>% mutate(medicaid_status = ifelse(medicaid_expansion == 1, "Expanded", "Not Expanded")) %>% dplyr::select(state_name, n_counties, mean_mortality, sd_mortality, medicaid_status) %>% kable(digits = 1, caption = "Average Mortality Rate by State") %>% kable_styling(bootstrap_options = c("striped", "hover")) %>% row_spec(which(state_means$medicaid_expansion == 1), background = "#E8F5E8")# Visualize the clusteringggplot(county_data, aes(x = reorder(state_name, mortality_rate), y = mortality_rate)) + geom_jitter(aes(color = factor(medicaid_expansion)), width = 0.2, alpha = 0.6, size = 2) + stat_summary(fun = mean, geom = "point", size = 4, shape = 23, fill = "red", color = "black") + scale_color_manual(values = c("0" = "coral", "1" = "steelblue"), labels = c("No Medicaid Expansion", "Medicaid Expanded"), name = "Policy Status") + labs(title = "County Mortality Rates Show Clear State Clustering", subtitle = "Red diamonds = state averages. Notice how counties cluster by state!", x = "State (ordered by mortality rate)", y = "Mortality Rate per 100,000") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "bottom")```**Figure Interpretation:** This visualization demonstrates the strong clustering pattern that necessitates multilevel modeling. Several key patterns emerge:1. **Clear state-level clustering**: Counties within each state cluster tightly around their state mean (red diamonds), with relatively small within-state variation compared to the large differences between states. This tight clustering within states violates the independence assumption of standard regression.2. **Systematic state differences**: States show dramatically different mortality patterns, with some states consistently experiencing mortality rates 100-200 deaths per 100,000 higher than others. The ordering from lowest to highest mortality reveals persistent state-level factors that affect all counties within a state's borders.3. **Policy pattern evidence**: The color coding reveals that Medicaid expansion states (blue) tend to cluster toward the lower end of the mortality distribution, while non-expansion states (coral) are more heavily represented among higher-mortality states. This suggests state-level policy decisions create systematic health advantages or disadvantages.This clustering pattern demonstrates why standard regression would be inappropriate—counties within the same state are clearly not independent observations, and ignoring this structure would lead to artificially precise estimates and incorrect conclusions about the significance of predictors.#### 2.2. Measuring Clustering: The ICCThe **Intracluster Correlation Coefficient (ICC)** tells us how much clustering exists. It answers: "How much of the total variance is between states vs. within states?"**Formula:** $ICC = \frac{\sigma^2_{between}}{\sigma^2_{between} + \sigma^2_{within}}$ Where:- $\sigma^2_{between}$ = how much states differ from each other- $\sigma^2_{within}$ = how much counties vary within each state#### 2.3. Calculating ICCThe easiest way is to fit an "empty" multilevel model with no predictors:```{r calculate-icc}# Step 1: Fit empty model (intercept + random state effects only)empty_model <- lmer(mortality_rate ~ 1 + (1|state_id), data = county_data)# Step 2: Look at the variance componentsprint(VarCorr(empty_model), comp = "Variance")# Step 3: Calculate ICC manuallyvariance_components <- as.data.frame(VarCorr(empty_model))between_state_var <- variance_components$vcov[1] # State variancewithin_state_var <- variance_components$vcov[2] # Residual variancetotal_variance <- between_state_var + within_state_varicc <- between_state_var / total_variance# Option 1: Clean table formatvariance_table <- data.frame( Component = c("Between-state variance", "Within-state variance", "Total variance"), Value = c(between_state_var, within_state_var, total_variance), Percentage = c(100*icc, 100*(1-icc), 100))kable(variance_table, digits = 1, col.names = c("Variance Component", "Value", "Percentage (%)"), caption = "Variance Decomposition Results") %>% kable_styling(bootstrap_options = c("striped", "hover")) %>% row_spec(3, bold = TRUE) %>% add_header_above(c(" " = 1, "Variance" = 2))cat(sprintf("\nIntraclass Correlation Coefficient (ICC) = %.3f (%.1f%%)", icc, 100*icc))```#### 2.4. Interpreting ICC Values**ICC Interpretation Guide:**- **ICC < 0.05**: Minimal clustering (multilevel modeling optional)- **ICC 0.05-0.10**: Small clustering (multilevel modeling recommended) - **ICC 0.10-0.25**: Moderate clustering (multilevel modeling important)- **ICC > 0.25**: Large clustering (multilevel modeling essential)**What ICC Means:**Our calculated ICC reveals striking insights about the structure of our health data. With an ICC of 0.798, this means:- 79.8% of mortality variance occurs BETWEEN states- 20.2% of mortality variance occurs WITHIN states - Counties in the same state are 79.8% more similar to each other than to random counties from other states**Decision for Our Analysis:**Since our ICC > 0.25 (and dramatically so), multilevel modeling is absolutely essential. This extremely high clustering indicates that ignoring the nested structure would lead to:- Severely underestimated standard errors (false precision)- Highly misleading conclusions about predictor significance- Complete failure to capture the dominant state-level effects on health outcomes**What This Tells Us About Health:**This ICC level reveals that state-level factors overwhelmingly dominate county health outcomes. Nearly four-fifths of the variation in mortality rates is attributable to systematic differences between states rather than local county characteristics. This suggests that state policies, healthcare systems, environmental regulations, and regional culture create extremely powerful contextual effects that strongly determine health outcomes for all counties within a state's borders.The relatively small within-state variation (20.2%) indicates that while county-level factors still matter, state context is by far the primary driver of health disparities. This multilevel structure demonstrates that effective health interventions must prioritize state-level policy changes, as local county-level interventions alone may have limited impact within the broader state policy environment.#### 2.5. What This Means for Our AnalysisThe ICC reveals important insights about our health data:1. **Substantial clustering exists**: Counties within states are much more similar than counties across states2. **State-level factors matter**: A significant portion of health outcomes is determined by state-level characteristics3. **Policy implications**: State policies (like Medicaid expansion) create systematic differences across all counties in a state4. **Methodological necessity**: We need multilevel models to properly account for this clusteringIn the next section, we'll build multilevel models that properly handle this nested structure and give us correct standard errors and meaningful policy insights.### 3. The Multilevel Model FrameworkMultilevel models (also called hierarchical linear models or mixed-effects models) are designed for **nested data structures**, where lower-level units (e.g., counties) are grouped within higher-level units (e.g., states). Standard regression assumes all observations are independent, but multilevel models account for **within-group correlation** and allow **group-level variation** in both intercepts and slopes.**Key concepts:**- **Levels:** - **Level 1 (County, $i$):** Individual observations at the lower level (counties within states). The **level 1 equation** models county-level mortality as: $$Mortality_{ij} = \beta_{0j} + \beta_1 Income_{ij} + \beta_2 Smoking_{ij} + \beta_3 Obesity_{ij} + \beta_4 Unemployment_{ij} + \cdots + r_{ij}$$ - **Level 2 (State, $j$):** Higher-level grouping units (states). The **level 2 equations** model how parameters vary across states: $$\beta_{0j} = \gamma_{00} + \gamma_{01} MedicaidExp_j + \gamma_{02} HealthSpend_j + u_{0j}$$ $$\beta_{1j} = \gamma_{10} + u_{1j} \quad \text{(if income effect varies by state)}$$- **Fixed effects** ($\gamma$ parameters): Average associations across all states, representing population-level relationships such as the mean effect of county income, smoking, or obesity on mortality.- **Random effects** ($u$ terms): State-specific deviations from the fixed effects, capturing heterogeneity across states: - $u_{0j} \sim N(0, \tau_{00})$: Random intercept capturing baseline mortality differences between states - $u_{1j} \sim N(0, \tau_{11})$: Random slope capturing how predictor effects vary between states- **Residuals** ($r_{ij} \sim N(0, \sigma^2)$): County-level deviations capturing **within-state variation** not explained by predictors.- **Cross-level interactions:** Terms where state-level factors modify county-level relationships: $$\text{Effect of Income} = \gamma_{10} + \gamma_{11} \times MedicaidExp_j$$ This captures whether Medicaid expansion strengthens or weakens the income-mortality association.- **Intraclass Correlation (ICC):** Proportion of total variance attributable to state-level clustering: $$ICC = \frac{\tau_{00}}{\tau_{00} + \sigma^2}$$ Higher ICC values indicate greater similarity within states relative to between-state differences.**Why use multilevel models?**1. **Corrects standard errors** for clustered data, providing appropriate uncertainty estimates2. **Partitions variation** into within-group and between-group components3. **Models contextual effects** by incorporating group-level predictors and cross-level interactions4. **Improves predictions** through partial pooling, borrowing strength across similar groups5. **Addresses ecological fallacy** by properly modeling relationships at multiple levelsThis framework enables us to examine how county characteristics relate to mortality while accounting for state-level clustering and investigating how state policies modify these relationships.### 4. Model Types and Complexity#### 4.1 Random Intercept Model (Simplest)The random intercept model allows only the **intercept** to vary by state, assuming that the effects of county-level predictors (income, smoking, obesity, etc.) are the same across all states.**Level 1 (County) Model:**$$Mortality_{ij} = \beta_{0j} + \beta_1 Income_{ij} + \beta_2 Smoking_{ij} + \beta_3 Obesity_{ij} + \beta_4 Unemployment_{ij} + \cdots + r_{ij}$$**Level 2 (State) Model:**$$\beta_{0j} = \gamma_{00} + u_{0j}$$$$\beta_1, \beta_2, \dots \text{ are fixed across states}$$**Combined (Mixed) Form:**$$Mortality_{ij} = \gamma_{00} + \beta_1 Income_{ij} + \beta_2 Smoking_{ij} + \beta_3 Obesity_{ij} + \beta_4 Unemployment_{ij} + \cdots + u_{0j} + r_{ij}$$**Variance Components:**- $u_{0j} \sim N(0, \tau_{00})$: Random intercept for state $j$, capturing deviations from the overall mean baseline mortality- $r_{ij} \sim N(0, \sigma^2)$: Residual for county $i$ in state $j$, capturing **within-state variation**- **Intraclass Correlation (ICC):**$$ICC = \frac{\tau_{00}}{\tau_{00} + \sigma^2}$$This ICC represents the proportion of total variance in mortality that occurs **between states** rather than within states, indicating the degree of clustering at the state level.```{r random-intercept-model}# Random intercept modelmodel_ri <- lmer(mortality_rate ~ median_income + smoking_rate + obesity_rate + unemployment_rate + hospitals_per_100k + primary_care_physicians_per_100k + factor(rural_status) + (1 | state_id), data = county_data, REML = TRUE)# Display resultstab_model(model_ri)# Calculate and display ICCicc_model <- performance::icc(model_ri)```The random intercept model shows:- **Fixed effects**: The average associations between predictors and mortality, pooled across all states. These correspond to the $\beta$ or $\gamma$ parameters in the model equations.- **Random intercept variance**: $\tau_{00} = 13387.26$ — the variation in baseline mortality **between states** captured by $u_{0j}$.- **Residual variance**: $\sigma^2 = 821.22$ — the remaining variation **within states** (between counties), represented by $r_{ij}$.- **Intraclass Correlation (ICC)**: 0.94 — the proportion of total variance in mortality attributable to differences **between states** after accounting for predictors.This high ICC value (0.94) indicates that even after accounting for county-level characteristics, the vast majority of remaining mortality variation occurs **between states** rather than within states. This suggests that unmeasured state-level factors (policies, healthcare systems, environmental conditions) dominate health outcomes, creating systematic advantages or disadvantages for all counties within a state's borders. County-level interventions alone may have limited impact without addressing these broader state-level determinants.```{r}# Visualize random intercept model# Create predictions with fixed slopes but varying interceptsstate_predictions_ri <- county_data %>% dplyr::select(state_id, state_name, median_income, mortality_rate) %>%group_by(state_id, state_name) %>%do(data.frame(median_income =seq(min(.$median_income), max(.$median_income), length.out =50),# Add the missing variables with their mean valuessmoking_rate =mean(county_data$smoking_rate),obesity_rate =mean(county_data$obesity_rate),unemployment_rate =mean(county_data$unemployment_rate),hospitals_per_100k =mean(county_data$hospitals_per_100k),primary_care_physicians_per_100k =mean(county_data$primary_care_physicians_per_100k),rural_status ="Suburban"# Use the most common category )) %>%ungroup()# Add predictions from random intercept modelstate_predictions_ri$predicted <-predict(model_ri, newdata = state_predictions_ri, allow.new.levels =FALSE)# Plot state-specific intercepts with parallel slopesp_intercepts <-ggplot() +# Individual county data pointsgeom_point(data = county_data, aes(x = median_income, y = mortality_rate, color = state_name),alpha =0.5, size =1.5) +# State-specific regression lines (parallel slopes)geom_line(data = state_predictions_ri, aes(x = median_income, y = predicted, color = state_name),linewidth =1.2) +labs(title ="State-Specific Baseline Mortality with Common Income Effect",subtitle ="Random intercept model: parallel lines show same income effect across states",x ="Median Household Income ($1000s)",y ="Predicted Mortality Rate (per 100,000)",color ="State") +theme_minimal() +theme(legend.position ="right")# Extract random interceptsrandom_intercepts <-ranef(model_ri)$state_idrandom_intercepts$state_name <- state_info$state_name[match(rownames(random_intercepts), state_info$state_id)]# Display random intercepts tablerandom_intercepts %>%arrange(`(Intercept)`) %>%mutate(`Intercept Deviation`=round(`(Intercept)`, 1) ) %>% dplyr::select(state_name, `Intercept Deviation`) %>%kable(caption ="State-Specific Random Intercepts") %>%kable_styling(bootstrap_options =c("striped", "hover", "condensed")) %>%column_spec(1, bold =TRUE)p_intercepts```**Figure Interpretation:** The random intercept model reveals important baseline differences across states while assuming uniform effects of county-level predictors. The parallel lines demonstrate several key patterns:1. **Intercept variation**: States have substantially different baseline mortality levels, with some states consistently higher or lower than others across all income levels. The vertical separation between lines represents these state-specific deviations from the national average.2. **Parallel slopes**: All lines have the same slope, indicating that the protective effect of income on mortality is assumed to be identical across all states. This constraint means that a $10,000 increase in median household income has the same predicted reduction in mortality rate regardless of the state context.3. **State rankings consistency**: States maintain their relative ranking across the income distribution. For example, if State A has higher baseline mortality than State B, this relationship persists at all income levels.The table of random intercepts quantifies these state-specific baseline deviations. Positive values indicate states with higher-than-average mortality rates after accounting for county-level characteristics, while negative values show states with lower-than-average mortality. These persistent state differences suggest the presence of unmeasured state-level factors (policies, healthcare systems, environmental conditions) that create systematic advantages or disadvantages for all counties within a state's borders.This pattern supports the need for state-level interventions to address the fundamental contextual factors driving these baseline differences, as county-level interventions alone cannot overcome the broader state policy environment.#### 4.2 Random Slope Model (More Complex)The random slope model allows the effect of median income to vary across states, recognizing that income may have different protective effects depending on state context:**Level 1 (County) Model:**$$Mortality_{ij} = \beta_{0j} + \beta_{1j} Income_{ij} + \beta_2 Smoking_{ij} + \beta_3 Obesity_{ij} + \beta_4 Unemployment_{ij} + \cdots + r_{ij}$$**Level 2 (State) Model:**$$\beta_{0j} = \gamma_{00} + u_{0j}$$$$\beta_{1j} = \gamma_{10} + u_{1j}$$$$\beta_2, \beta_3, \dots \text{ remain fixed across states}$$**Combined (Mixed) Form:**$$Mortality_{ij} = \gamma_{00} + \gamma_{10} Income_{ij} + \beta_2 Smoking_{ij} + \beta_3 Obesity_{ij} + \cdots + u_{0j} + u_{1j} Income_{ij} + r_{ij}$$Where:- $u_{0j}$: Random intercept for state $j$ (baseline mortality differences between states)- $u_{1j}$: Random slope for median income in state $j$ (how the income-mortality relationship varies by state)This extension allows us to test whether the protective effect of higher income is consistent across all states or varies systematically based on state-level characteristics.```{r random-slope-model}#| warning: false# Random slope model - allowing income effect to vary by statemodel_rs <- lmer(mortality_rate ~ median_income + smoking_rate + obesity_rate + unemployment_rate + hospitals_per_100k + primary_care_physicians_per_100k + factor(rural_status) + (1 + median_income | state_id), data = county_data, REML = TRUE)# Display resultstab_model(model_rs)# Extract random effectsrandom_effects <- ranef(model_rs)$state_idrandom_effects$state_name <- state_info$state_name[match(rownames(random_effects), state_info$state_id)]```The random slope model shows:- **Fixed effects**: The average relationships between predictors and mortality across all states.- **Random intercept variance**: $\tau_{00} = 4529.16$ — variation in baseline mortality **between states** (reduced from 13387.26 in random intercept model).- **Random slope variance (median income)**: $\tau_{11} = 0.98$ — variation in the effect of median income **across states**.- **Intercept–slope correlation**: $\rho_{01} = 0.64$ — positive correlation indicating that states with higher baseline mortality tend to have stronger (more negative) income effects.- **Residual variance**: $\sigma^2 = 651.47$ — the remaining within-state variation (between counties).The model shows improved fit over the random intercept model (Conditional R² increased from 0.951 to 0.962). **Variance interpretation**: The substantial reduction in random intercept variance (from 13387 to 4529) occurs because allowing income effects to vary by state captures some of the between-state differences that were previously attributed to baseline differences. States with stronger income-mortality relationships now have this heterogeneity explicitly modeled rather than absorbed into their random intercepts. The positive correlation (0.64) suggests that in states where mortality is generally higher, income differences may be even more consequential for health outcomes.```{r visualize-random-slopes, fig.width=12, fig.height=8}# Create state-specific regression linesstate_predictions <- county_data %>% dplyr::select(state_id, state_name, median_income, mortality_rate) %>% group_by(state_id, state_name) %>% do(data.frame( median_income = seq(min(.$median_income), max(.$median_income), length.out = 50), # Add the missing variables with their mean values smoking_rate = mean(county_data$smoking_rate), obesity_rate = mean(county_data$obesity_rate), unemployment_rate = mean(county_data$unemployment_rate), hospitals_per_100k = mean(county_data$hospitals_per_100k), primary_care_physicians_per_100k = mean(county_data$primary_care_physicians_per_100k), rural_status = "Suburban" # Use the most common category )) %>% ungroup()# Add predictions from random slope modelstate_predictions$predicted <- predict(model_rs, newdata = state_predictions, allow.new.levels = FALSE)# Plot state-specific slopesp_slopes <- ggplot() + # Individual county data points geom_point(data = county_data, aes(x = median_income, y = mortality_rate, color = state_name), alpha = 0.5, size = 1.5) + # State-specific regression lines geom_line(data = state_predictions, aes(x = median_income, y = predicted, color = state_name), linewidth = 1.2) + labs(title = "State-Specific Income-Mortality Relationships", subtitle = "Random slope model allows income effects to vary by state", x = "Median Household Income ($1000s)", y = "Predicted Mortality Rate (per 100,000)", color = "State") + theme_minimal() + theme(legend.position = "right")# Display random effects tablerandom_effects %>% arrange(median_income) %>% mutate( `Intercept Deviation` = round(`(Intercept)`, 1), `Income Slope Deviation` = round(median_income, 3) ) %>% dplyr::select(state_name, `Intercept Deviation`, `Income Slope Deviation`) %>% kable(caption = "State-Specific Random Effects") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% column_spec(1, bold = TRUE)p_slopes```**Figure Interpretation:** The random slope model reveals important heterogeneity across states. Each colored line represents a different state's income-mortality relationship. We can see that:1. **Slope variation**: Some states (e.g., those with steeper negative slopes) show stronger protective effects of income on mortality, while others show weaker relationships.2. **Intercept variation**: States start at different baseline mortality levels (where lines would intersect the y-axis).3. **Policy implications**: This heterogeneity suggests that income-support policies might have different effectiveness across states, and one-size-fits-all federal policies may not be optimal.The table of random effects quantifies these state-specific deviations from the average relationship, helping identify which states might benefit most from targeted interventions.#### 4.3 Cross-Level Interaction Model (Most Complex)The cross-level interaction model incorporates state-level predictors and examines how state policies modify the effects of county-level variables. We introduce this complexity gradually, first adding state-level predictors, then their interactions with county-level variables.**State-Level Predictors Added:**- Medicaid expansion status (policy variable)- State health spending per capita - Tobacco tax rates**Level 1 (County) Model:**$$Mortality_{ij} = \beta_{0j} + \beta_{1j} Income_{ij} + \beta_{2j} Smoking_{ij} + \beta_3 Obesity_{ij} + \beta_4 Unemployment_{ij} + \cdots + r_{ij}$$**Level 2 (State) Model:**$$\beta_{0j} = \gamma_{00} + \gamma_{01} MedicaidExp_j + \gamma_{02} HealthSpend_j + \gamma_{03} TobaccoTax_j + u_{0j}$$$$\beta_{1j} = \gamma_{10} + \gamma_{11} MedicaidExp_j$$$$\beta_{2j} = \gamma_{20} + \gamma_{21} TobaccoTax_j$$$$\beta_3, \beta_4, \dots \text{ remain fixed across states}$$**Combined (Mixed) Form:**$$\begin{align}Mortality_{ij} = &\gamma_{00} + \gamma_{10} Income_{ij} + \gamma_{20} Smoking_{ij} + \gamma_{01} MedicaidExp_j + \gamma_{02} HealthSpend_j + \gamma_{03} TobaccoTax_j \\&+ \gamma_{11} (Income_{ij} \times MedicaidExp_j) + \gamma_{21} (Smoking_{ij} \times TobaccoTax_j) + \cdots + u_{0j} + r_{ij}\end{align}$$This model captures:- **Main effects** of both county-level and state-level predictors- **Cross-level interactions** showing how state policies modify county-level relationships- $u_{0j}$: Random intercept capturing remaining baseline variation across states after accounting for state-level predictors```{r cross-level-model}# Cross-level interaction modelmodel_cli <- lmer(mortality_rate ~ # County-level main effects median_income + smoking_rate + obesity_rate + unemployment_rate + hospitals_per_100k + primary_care_physicians_per_100k + factor(rural_status) + # State-level main effects medicaid_expansion + state_health_spending + tobacco_tax + # Cross-level interactions median_income:medicaid_expansion + smoking_rate:tobacco_tax + # Random effects (1 | state_id), data = county_data, REML = TRUE)tab_model(model_cli)```The cross-level interaction model shows:- **Fixed effects**: - Average county-level associations between predictors (e.g., income, smoking, obesity) and mortality. - State-level effects (e.g., Medicaid expansion, state health spending, tobacco tax) on mortality. - **Cross-level interactions**: how state policies modify county-level relationships (e.g., whether income has a stronger protective effect in Medicaid expansion states, or smoking is more harmful in states with low tobacco taxes).- **Random intercept variance**: $\tau_{00} = 287.93$ — variation in baseline mortality **between states** after including state-level predictors.- **Residual variance**: $\sigma^2 = 603.88$ — remaining variation **within states** (between counties).- **Intraclass Correlation (ICC)**: 0.32 — the proportion of total variance in mortality that lies **between states** after accounting for both county- and state-level predictors.**Variance interpretation**: The dramatic reduction in ICC from 0.94 (random intercept) to 0.32 (cross-level interaction) occurs because state-level predictors (Medicaid expansion, health spending, tobacco taxes) explain much of the between-state variation that was previously unexplained. The inclusion of these measured state characteristics reduces the remaining unobserved state heterogeneity from 94% to 32% of total variance. This suggests that state policies and their interactions with local characteristics are primary drivers of health disparities across the United States. The cross-level interactions further explain why certain county-level factors (like income) have different effects in different policy contexts.```{r interaction-visualization, fig.width=12, fig.height=8}# Visualize cross-level interactions# 1. Income effect by Medicaid expansion statusincome_range <- seq(min(county_data$median_income), max(county_data$median_income), length.out = 50)pred_data <- expand.grid( median_income = income_range, medicaid_expansion = c(0, 1), smoking_rate = mean(county_data$smoking_rate), obesity_rate = mean(county_data$obesity_rate), unemployment_rate = mean(county_data$unemployment_rate), hospitals_per_100k = mean(county_data$hospitals_per_100k), primary_care_physicians_per_100k = mean(county_data$primary_care_physicians_per_100k), rural_status = "Suburban", state_health_spending = mean(county_data$state_health_spending), tobacco_tax = mean(county_data$tobacco_tax), state_id = 1 # Use a reference state for prediction)pred_data$predicted <- predict(model_cli, newdata = pred_data, re.form = NA)p_interaction1 <- ggplot(pred_data, aes(x = median_income, y = predicted, color = factor(medicaid_expansion))) + geom_line(linewidth = 1.5) + scale_color_manual(values = c("0" = "red", "1" = "blue"), labels = c("No Medicaid Expansion", "Medicaid Expanded"), name = "Policy Status") + labs(title = "Cross-Level Interaction: Income × Medicaid Expansion", subtitle = "Medicaid expansion modifies the income-mortality relationship", x = "Median Household Income ($1000s)", y = "Predicted Mortality Rate (per 100,000)") + theme_minimal() + theme(legend.position = "bottom")# 2. Smoking effect by tobacco tax levelssmoking_range <- seq(min(county_data$smoking_rate), max(county_data$smoking_rate), length.out = 50)tobacco_levels <- quantile(county_data$tobacco_tax, c(0.25, 0.75))pred_data2 <- expand.grid( smoking_rate = smoking_range, tobacco_tax = tobacco_levels, median_income = mean(county_data$median_income), medicaid_expansion = 0, obesity_rate = mean(county_data$obesity_rate), unemployment_rate = mean(county_data$unemployment_rate), hospitals_per_100k = mean(county_data$hospitals_per_100k), primary_care_physicians_per_100k = mean(county_data$primary_care_physicians_per_100k), rural_status = "Suburban", state_health_spending = mean(county_data$state_health_spending), state_id = 1)pred_data2$predicted <- predict(model_cli, newdata = pred_data2, re.form = NA)p_interaction2 <- ggplot(pred_data2, aes(x = smoking_rate, y = predicted, color = factor(tobacco_tax))) + geom_line(linewidth = 1.5) + scale_color_manual(values = c("steelblue", "orange"), labels = c("Low Tobacco Tax", "High Tobacco Tax"), name = "Tax Policy") + labs(title = "Cross-Level Interaction: Smoking Rate × Tobacco Tax", subtitle = "State tobacco policies modify the smoking-mortality relationship", x = "County Smoking Rate (%)", y = "Predicted Mortality Rate (per 100,000)") + theme_minimal() + theme(legend.position = "bottom")p_interaction1p_interaction2```**Figure Interpretation:** The cross-level interaction results reveal important policy insights:**Top Panel (Income × Medicaid Expansion):** The diverging lines show that the protective effect of income on mortality differs by Medicaid expansion status. In states without Medicaid expansion (red line), the income gradient is steeper—meaning low-income counties have particularly high mortality. In expansion states (blue line), this gradient is flatter, suggesting that Medicaid expansion partially compensates for income-based health disparities. This finding supports the hypothesis that Medicaid expansion provides a safety net that reduces health inequalities.**Bottom Panel (Smoking × Tobacco Tax):** Higher tobacco taxes (orange line) reduce the harmful effect of smoking on mortality compared to low-tax states (blue line). The flatter slope in high-tax states suggests that tobacco taxation not only reduces smoking rates but also mitigates the health consequences of smoking, possibly through reduced intensity of smoking or funding for cessation programs. This demonstrates how state policies can modify behavioral risk factors' impacts on population health.#### 4.4 Model Selection Decision FrameworkChoosing the appropriate level of model complexity requires balancing theoretical justification, statistical considerations, and practical constraints. The literature provides varied guidance on sample size requirements, reflecting the complexity of multilevel modeling decisions.**Step 1: Start Simple - Random Intercept Model**- **Always begin here** regardless of your research question- Provides baseline ICC and model fit statistics- Establishes whether multilevel modeling is necessary (ICC > 0.05)**Step 2: Assess Need for Random Slopes***Primary Consideration: Theoretical Justification*- Do you have strong theory suggesting predictor effects vary meaningfully across groups?- Are cross-group differences in slopes substantively important for your research question?- Do existing studies suggest heterogeneous effects in your domain?*Statistical Evidence*- Is there significant variability in the predictor across Level 2 units?- Do exploratory plots suggest non-parallel relationships?- Does likelihood ratio test support added complexity?*Sample Size Considerations*The literature provides conflicting guidance on minimum group numbers:- **Conservative estimates**: 30+ groups for reliable variance component estimation- **Moderate estimates**: 10-20 groups may be sufficient with large within-group samples- **Liberal estimates**: As few as 5-8 groups in some simulation studies**Reality**: Your decision should prioritize theoretical justification over arbitrary sample size rules.**Step 3: Consider Cross-Level Interactions***Theoretical Foundation (Most Important)*- Do you have measured Level 2 variables that theoretically moderate Level 1 relationships?- Are policy or contextual effects central to your research question?- Is there substantive reason to believe effects vary systematically (not just randomly)?*Model Building Strategy*- Add Level 2 main effects before interactions- Include interactions only with strong theoretical rationale- Consider whether random slopes are needed before cross-level interactions**Decision Framework Table:**```{r model-selection-table}# Create decision framework tabledecision_framework <- data.frame( `Model Type` = c("Random Intercept", "Random Slope", "Cross-Level Interaction"), `When to Use` = c( "Default starting point; Interest in average effects; Theory assumes constant effects", "Theory suggests effect heterogeneity; Exploratory evidence of variation; Policy questions about differential impact", "Specific hypotheses about policy/context moderation; Need to explain slope variation" ), `Sample Considerations` = c( "Flexible - works with small Level 2 samples", "Literature varies: 10-30+ groups depending on context and within-group size", "Similar to random slopes, plus need variation in Level 2 moderators" ), `Key Trade-offs` = c( "Simple but may miss important heterogeneity; Assumes one-size-fits-all", "More realistic but convergence issues; Harder to interpret", "Policy-relevant but complex; Risk of overfitting and interpretation challenges" ), `Health Research Example` = c( "Average effect of smoking on mortality (assumes same across all states)", "Smoking effects vary by state (captures heterogeneity without explaining why)", "Tobacco tax policy moderates smoking-mortality relationship (explains why effects vary)" ))decision_framework %>% kable(caption = "Model Selection Decision Framework") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% column_spec(1, bold = TRUE, width = "15%") %>% column_spec(2, width = "30%") %>% column_spec(3, width = "25%") %>% column_spec(4, width = "30%")```**Practical Guidelines:***What the Literature Actually Shows*- Random intercept models: Generally robust even with small Level 2 samples- Random slope models: Performance depends on Level 2 sample size, within-group sample sizes, and true effect heterogeneity- Cross-level interactions: Require sufficient variation in both Level 1 and Level 2 variables*Pragmatic Approach*1. **Theory first**: Strong theoretical justification trumps sample size concerns2. **Convergence as reality check**: Non-convergent models signal overly complex specifications3. **Substantive significance**: Focus on effect sizes and policy implications, not just statistical significance4. **Sensitivity analysis**: Test robustness of conclusions across different model specifications**Recommended Workflow:**1. **Fit random intercept model** → examine ICC and theoretical fit2. **Test random slopes** for theoretically important predictors → check convergence and improvement3. **Add Level 2 predictors** to explain group-level variation → assess explanatory power4. **Consider cross-level interactions** only with strong theoretical basis → focus on policy/practical relevance5. **Model comparison**: Use multiple criteria (theory, fit, convergence, interpretability)6. **Sensitivity analysis**: Test key conclusions across different specifications### 5. Model Comparison, Diagnostics, and Best Practices#### 5.1 Model Comparison```{r model-comparison}# Compare models using information criteriamodel_comparison <- data.frame( Model = c("Random Intercept", "Random Slope", "Cross-Level Interaction"), AIC = c(AIC(model_ri), AIC(model_rs), AIC(model_cli)), BIC = c(BIC(model_ri), BIC(model_rs), BIC(model_cli)), `Log Likelihood` = c(logLik(model_ri), logLik(model_rs), logLik(model_cli)), `Parameters` = c(attr(logLik(model_ri), "df"), attr(logLik(model_rs), "df"), attr(logLik(model_cli), "df")))model_comparison %>% mutate( AIC = round(AIC, 1), BIC = round(BIC, 1), Log.Likelihood = round(Log.Likelihood, 1) ) %>% kable(caption = "Model Comparison for Health Data") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% column_spec(1, bold = TRUE)# Likelihood ratio testsanova(model_ri, model_rs, model_cli)```The model comparison shows that increasing complexity (from random intercept to random slope to cross-level interactions) improves model fit as indicated by lower AIC values. The likelihood ratio tests confirm that each additional complexity is statistically justified.#### 5.2 Model Diagnostics```{r diagnostics, fig.width=12, fig.height=10}# 1. Residual diagnosticscounty_data$residuals <- residuals(model_cli)county_data$fitted <- fitted(model_cli)# 2. Random effects diagnosticsrandom_intercepts <- ranef(model_cli)$state_id[,1]# Create diagnostic plotsp_residuals <- ggplot(county_data, aes(x = fitted, y = residuals)) + geom_point(alpha = 0.5) + geom_hline(yintercept = 0, color = "red", linetype = "dashed") + geom_smooth(method = "loess", se = FALSE, color = "blue") + labs(title = "Residuals vs Fitted Values", subtitle = "Check for homoscedasticity and linearity", x = "Fitted Values", y = "Residuals") + theme_minimal()p_qq <- ggplot(county_data, aes(sample = residuals)) + stat_qq() + stat_qq_line(color = "red") + labs(title = "Q-Q Plot of Residuals", subtitle = "Check normality assumption", x = "Theoretical Quantiles", y = "Sample Quantiles") + theme_minimal()p_random <- data.frame(random_intercepts = random_intercepts) %>% ggplot(aes(sample = random_intercepts)) + stat_qq() + stat_qq_line(color = "red") + labs(title = "Q-Q Plot of Random Intercepts", subtitle = "Check normality of random effects", x = "Theoretical Quantiles", y = "Random Intercepts") + theme_minimal()p_leverage <- ggplot(county_data, aes(x = fitted, y = sqrt(abs(residuals)))) + geom_point(alpha = 0.5) + geom_smooth(method = "loess", se = FALSE, color = "blue") + labs(title = "Scale-Location Plot", subtitle = "Check homoscedasticity", x = "Fitted Values", y = "√|Residuals|") + theme_minimal()p_residualsp_qqp_randomp_leverage```**Figure Interpretation:** The diagnostic plots assess key model assumptions:1. **Residuals vs Fitted (top left)**: The blue loess line should be relatively flat around zero. Minor deviations suggest the model captures the main patterns well, though some non-linearity may remain.2. **Q-Q Plot of Residuals (top right)**: Points following the red line indicate normally distributed residuals. Deviations at the tails are common in health data but shouldn't invalidate main conclusions.3. **Q-Q Plot of Random Intercepts (bottom left)**: Checks if state-level random effects are normally distributed. With only 8 states, perfect normality is unlikely, but major deviations would be concerning.4. **Scale-Location Plot (bottom right)**: A horizontal trend would indicate constant variance. Some heteroscedasticity is visible but not severe enough to require transformation.#### 5.3 Variance Decomposition Analysis```{r variance-decomposition, fig.width=12, fig.height=6}# Extract variance components from different modelsvar_components <- data.frame( Model = c("Null Model", "Random Intercept", "Random Slope", "Cross-Level"), State_Variance = c( VarCorr(lmer(mortality_rate ~ 1 + (1|state_id), county_data))$state_id[1,1], VarCorr(model_ri)$state_id[1,1], VarCorr(model_rs)$state_id[1,1], VarCorr(model_cli)$state_id[1,1] ), County_Variance = c( attr(VarCorr(lmer(mortality_rate ~ 1 + (1|state_id), county_data)), "sc")^2, attr(VarCorr(model_ri), "sc")^2, attr(VarCorr(model_rs), "sc")^2, attr(VarCorr(model_cli), "sc")^2 ))var_components <- var_components %>% mutate( Total_Variance = State_Variance + County_Variance, Prop_State = State_Variance / Total_Variance, Prop_County = County_Variance / Total_Variance )# Create visualizationlibrary(tidyr)var_long <- var_components %>% dplyr::select(Model, Prop_State, Prop_County) %>% pivot_longer(cols = c(Prop_State, Prop_County), names_to = "Level", values_to = "Proportion") %>% mutate( Level = gsub("Prop_", "", Level), Model = factor(Model, levels = c("Null Model", "Random Intercept", "Random Slope", "Cross-Level")) )ggplot(var_long, aes(x = Model, y = Proportion, fill = Level)) + geom_col(position = "stack", alpha = 0.8) + scale_fill_manual(values = c("State" = "#E74C3C", "County" = "#3498DB")) + scale_y_continuous(labels = scales::percent) + labs(title = "Variance Decomposition Across Model Types", subtitle = "How much variation is explained at each level", x = "Model Type", y = "Proportion of Variance", fill = "Level") + theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1))```**Figure Interpretation:** This variance decomposition reveals the dramatic impact of including state-level predictors and interactions on our understanding of health disparities. The progression shows four key insights:1. **Null Model**: Nearly 80% of mortality variation occurs between states, with only 20% within states—demonstrating massive state-level clustering that demands multilevel analysis.2. **Random Intercept Model**: Adding county-level predictors (income, smoking, healthcare access) barely changes the variance partition. State-level factors still dominate, indicating that county characteristics alone cannot explain why some states consistently outperform others.3. **Random Slope Model**: Allowing income effects to vary by state maintains similar proportions, suggesting that heterogeneous county-level relationships don't substantially alter the fundamental state-county variance structure.4. **Cross-Level Interaction Model**: The dramatic reversal occurs when we include state policies and their interactions—now 70% of variance is within-state and only 30% between-state. This transformation reveals that measured state-level factors (Medicaid expansion, health spending, tobacco taxes) explain most of the systematic differences between states.### 6. Centering in Multilevel Models: A Critical DecisionOne of the most important but often overlooked decisions in multilevel modeling is how to center continuous predictors. The choice of centering strategy fundamentally affects the interpretation of coefficients and can lead to dramatically different conclusions about the same data.#### 6.1 The Centering ProblemWhen we introduce a continuous predictor in a multilevel model, we have to decide how to **center** it. This choice changes the meaning of the intercept and sometimes the slope. Here are the three main approaches, using median county income as an example:1. **Uncentered (raw values)** * We just use the variable as it is: $X_{ij} = Income_{ij}$ * **Interpretation**: The intercept represents the predicted outcome when income equals zero, which may be unrealistic (few counties have $0 income). * **Use case**: Sometimes helpful if the zero point is substantively meaningful (e.g., number of children = 0).2. **Grand mean centering** * Subtract the **overall mean** (across all counties) from each county's value: $X_{ij} = Income_{ij} - \overline{Income}$ * **Interpretation**: The intercept now represents the predicted outcome for a county with **average income**. Slopes represent the effect of income relative to the national average. * **Use case**: This is often the most interpretable default, because it avoids weird "zero" values and makes coefficients easier to explain across the whole sample.3. **Group mean centering (a.k.a. within-group centering)** * Subtract the **state mean** from each county's value: $X_{ij} = Income_{ij} - \overline{Income_j}$ * **Interpretation**: The intercept now represents the predicted outcome for a county at the **average income level of its own state**. The slope shows how counties differ from each other *within the same state*, independent of between-state differences. * **Use case**: This is crucial when we want to isolate **within-state effects** and avoid conflating them with between-state variation.Where $\overline{Income}$ is the grand mean across all counties and $\overline{Income_j}$ is the mean income within state $j$.```{r centering-demonstration}# Calculate different centering approachescounty_data <- county_data %>% mutate( # Original income (uncentered) income_raw = median_income, # Grand mean centering income_gmc = median_income - mean(median_income), # Group mean centering (within state) income_cwc = ave(median_income, state_id, FUN = function(x) x - mean(x)) ) %>% group_by(state_id) %>% mutate( # State mean income for comparison state_mean_income = mean(median_income) ) %>% ungroup()# Show the differencescentering_example <- county_data %>% filter(state_id %in% c(1, 2)) %>% dplyr::select(county_id, state_name, income_raw, income_gmc, income_cwc, state_mean_income) %>% slice_head(n = 6)centering_example %>% mutate_if(is.numeric, round, 2) %>% kable(caption = "Centering Approaches for Median Income", col.names = c("County", "State", "Raw Income", "Grand Mean Centered", "Group Mean Centered", "State Mean")) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed"))``````{r contextual-effects-setup}# Create between and within components for contextual effectscounty_data <- county_data %>% group_by(state_id) %>% mutate( # Between component: state mean (contextual effect) income_between = mean(median_income), # Within component: deviation from state mean income_within = median_income - mean(median_income), # Also create these for other key variables smoking_between = mean(smoking_rate), smoking_within = smoking_rate - mean(smoking_rate), obesity_between = mean(obesity_rate), obesity_within = obesity_rate - mean(obesity_rate) ) %>% ungroup()``````{r contextual-visualization, fig.width=14, fig.height=10}# Create comprehensive visualization of contextual effectsset.seed(789)# Panel 1: Show the decomposition visuallyp_decomp <- county_data %>% filter(state_id %in% c(1, 2, 4)) %>% ggplot(aes(x = median_income, y = mortality_rate)) + geom_point(aes(color = state_name), size = 2.5, alpha = 0.6) + geom_vline(aes(xintercept = income_between, color = state_name), linetype = "dashed", linewidth = 1) + geom_smooth(aes(group = state_name, color = state_name), method = "lm", se = FALSE, linewidth = 1) + geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 1.5, linetype = "dotted") + labs(title = "A. Decomposing Total, Between, and Within Effects", subtitle = "Vertical lines = state means; Colored lines = within-state; Black dotted = overall", x = "Median Income ($1000s)", y = "Mortality Rate (per 100,000)") + theme_minimal() + theme(legend.position = "bottom")# Panel 2: Between-state effects (contextual)state_means <- county_data %>% group_by(state_id, state_name) %>% summarize( mean_income = mean(median_income), mean_mortality = mean(mortality_rate), mean_smoking = mean(smoking_rate), medicaid = first(medicaid_expansion), .groups = "drop" )p_between <- ggplot(state_means, aes(x = mean_income, y = mean_mortality)) + geom_point(aes(size = 3, color = factor(medicaid)), alpha = 0.8) + geom_smooth(method = "lm", se = TRUE, color = "darkblue", linewidth = 1.2) + geom_text(aes(label = state_name), vjust = -1, hjust = 0.5, size = 3) + scale_color_manual(values = c("0" = "red", "1" = "blue"), labels = c("No Expansion", "Medicaid Expanded"), name = "Medicaid Status") + labs(title = "B. Between-State (Contextual) Effect", subtitle = "State-level relationship between mean income and mean mortality", x = "State Mean Income ($1000s)", y = "State Mean Mortality Rate") + theme_minimal() + theme(legend.position = "bottom") + guides(size = "none")# Panel 3: Within-state effects pooledp_within <- county_data %>% ggplot(aes(x = income_within, y = mortality_rate - ave(mortality_rate, state_id))) + geom_point(aes(color = state_name), alpha = 0.4, size = 1.5) + geom_smooth(method = "lm", se = TRUE, color = "black", linewidth = 1.5) + geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) + geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) + labs(title = "C. Within-State Effect (Pooled)", subtitle = "County deviations from state means", x = "Income - State Mean ($1000s)", y = "Mortality - State Mean") + theme_minimal() + theme(legend.position = "none")# Combine panelsp_decompp_betweenp_within```**Figure Interpretation:** This three-panel visualization highlights why it is essential to distinguish within- and between-state effects.* **Panel A** illustrates how the overall association (black dotted line) blends two sources of variation: differences **within states** (colored regression lines) and differences **between states** (vertical dashed lines marking state means). Each state not only has its own income–mortality slope, but states also vary in their average levels of income and mortality.* **Panel B** depicts the **between-state effect**: states with higher average income tend to experience lower average mortality. The contrast between red (no Medicaid expansion) and blue (expansion) states suggests that policy context shapes outcomes even at similar income levels. This is the ecological, or contextual, relationship operating at the state level.* **Panel C** isolates the **within-state effect** by centering counties on their state means. This shows how counties diverge from their state’s average income and mortality. The negative slope reveals that, within a given state, wealthier counties consistently experience lower mortality—a relationship that holds once state-level differences are removed.#### 6.2 Estimating and Interpreting Contextual Effects ModelsThe contextual effects model explicitly separates within-group and between-group effects, providing the most complete understanding of how predictors operate at different levels. This approach is the most sophisticated form of centering because it simultaneously includes **both** the group-mean centered variable (within-group effect) and the group mean itself (between-group effect) as separate predictors. This **dual-centering** strategy allows researchers to test whether the effect of income differs when comparing counties within the same state versus comparing states with different average income levels.The key question the contextual effects model addresses is: **"Does living in a high-income state provide health benefits beyond what we'd expect from individual county income alone?"**```{r contextual-models}# Model 1: Traditional random intercept model (for comparison)model_traditional <- lmer(mortality_rate ~ median_income + smoking_rate + obesity_rate + (1 | state_id), data = county_data, REML = TRUE)# Model 2: Full contextual effects modelmodel_contextual <- lmer(mortality_rate ~ # Within-state effects (group-mean centered) income_within + smoking_within + obesity_within + # Between-state effects (contextual/ecological) income_between + smoking_between + obesity_between + # Random intercept (1 | state_id), data = county_data, REML = TRUE)tab_model(model_traditional, model_contextual, title = "Comparing Traditional vs. Contextual Effects Models")```The comparison between traditional and contextual effects models reveals striking differences in how we understand income's relationship with mortality:**Traditional Model Findings:**- **Median income coefficient**: -2.55 (p < 0.001) - Each $1,000 increase in county income associates with 2.55 fewer deaths per 100,000- **Strong statistical significance** suggests income clearly matters for health outcomes- **High ICC (0.89)** indicates substantial clustering within states**Contextual Effects Model Findings:**- **Within-state income effect**: -2.55 (p < 0.001) - Identical to the traditional model, showing robust within-state relationships- **Between-state income effect**: -19.26 (p = 0.719) - Large magnitude but not statistically significant- **Wide confidence intervals** (-124.20 to 85.68) for between-state effects reflect uncertainty with only 8 states**Critical Interpretation Issues:**The contextual effects results suggest a potentially important pattern—the between-state coefficient is nearly 8 times larger than the within-state effect—but our analysis faces a fundamental limitation. With only 8 states in the dataset, we lack sufficient power to detect between-state relationships reliably. The wide confidence intervals and non-significant p-values for all between-state effects (income, smoking, obesity) likely reflect this small sample size rather than true absence of contextual effects.**What This Means for Health Research:**1. **Within-state relationships are robust**: The consistent -2.55 coefficient across models shows that county-level income differences within states reliably predict mortality differences.2. **Between-state effects are underpowered**: The large but non-significant between-state coefficients suggest potential contextual effects that cannot be reliably estimated with 8 states.3. **Methodological caution needed**: This demonstrates why contextual effects models require substantial numbers of higher-level units (typically 20+ states/regions) for meaningful between-group comparisons.The analysis illustrates both the promise and limitations of contextual effects modeling in health research—while theoretically powerful, the approach demands adequate sample sizes at all levels to produce reliable conclusions about how context modifies individual-level relationships.```{r contextual-test}# Extract coefficients for testingwithin_income <- fixef(model_contextual)["income_within"]between_income <- fixef(model_contextual)["income_between"]contextual_effect_income <- between_income - within_incomewithin_smoking <- fixef(model_contextual)["smoking_within"]between_smoking <- fixef(model_contextual)["smoking_between"]contextual_effect_smoking <- between_smoking - within_smoking# Create comprehensive results tablecontextual_results <- data.frame( Variable = c("Income", "Income", "Income", "Smoking", "Smoking", "Smoking"), Component = c("Within-State Effect", "Between-State Effect", "Contextual Effect (Difference)", "Within-State Effect", "Between-State Effect", "Contextual Effect (Difference)"), Coefficient = c( round(within_income, 3), round(between_income, 3), round(contextual_effect_income, 3), round(within_smoking, 3), round(between_smoking, 3), round(contextual_effect_smoking, 3) ), Interpretation = c( "Effect of $1K income increase within a state", "Effect of living in a state with $1K higher mean income", "Additional contextual benefit of high-income state environment", "Effect of 1% smoking increase within a state", "Effect of living in a state with 1% higher mean smoking", "Additional contextual harm of high-smoking state environment" ))contextual_results %>% kable(caption = "Contextual Effects Analysis: Do State Contexts Matter Beyond Individual Characteristics?") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% column_spec(1, bold = TRUE) %>% column_spec(2, bold = TRUE) %>% row_spec(c(3, 6), background = "#FFF3CD") %>% pack_rows("Income Effects", 1, 3, label_row_css = "background-color: #E8F4F8;") %>% pack_rows("Smoking Effects", 4, 6, label_row_css = "background-color: #F8E8E8;")cat(sprintf("\n**Key Findings:**"))cat(sprintf("\n• Income contextual effect: %.3f (%.1f%% stronger between-state effect)", contextual_effect_income, 100*(contextual_effect_income/abs(within_income))))cat(sprintf("\n• Smoking contextual effect: %.3f (%.1f%% difference in between-state effect)", contextual_effect_smoking, 100*(abs(contextual_effect_smoking)/abs(within_smoking))))```The critical test in contextual effects models is whether the **within-group** and **between-group** coefficients differ significantly. If they're the same, there's no contextual effect - state context doesn't matter beyond individual county characteristics.### 7. Practical Applications and Policy Implications#### 7.1 Policy Evaluation Framework```{r policy-evaluation}# Simulate policy intervention: What if all states expanded Medicaid?county_data_counterfactual <- county_data %>% mutate(medicaid_expansion_counterfactual = 1)# Predict mortality under universal Medicaid expansionpred_universal <- predict(model_cli, newdata = county_data_counterfactual, re.form = NULL)pred_current <- predict(model_cli, newdata = county_data, re.form = NULL)# Calculate policy impactpolicy_impact <- county_data %>% mutate( current_predicted = pred_current, universal_predicted = pred_universal, mortality_reduction = current_predicted - universal_predicted, lives_saved = (mortality_reduction / 100000) * population )# Summarize impact by statestate_impact <- policy_impact %>% filter(medicaid_expansion == 0) %>% # Only non-expansion states would be affected group_by(state_name) %>% summarize( counties_affected = n(), total_population = sum(population), avg_mortality_reduction = mean(mortality_reduction), total_lives_saved = sum(lives_saved), .groups = "drop" ) %>% arrange(desc(total_lives_saved))state_impact %>% mutate( total_population = round(total_population / 1000000, 2), avg_mortality_reduction = round(avg_mortality_reduction, 1), total_lives_saved = round(total_lives_saved, 0) ) %>% kable(caption = "Projected Impact of Universal Medicaid Expansion", col.names = c("State", "Counties Affected", "Population (Millions)", "Avg Mortality Reduction", "Lives Saved Annually")) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% column_spec(1, bold = TRUE)```This policy simulation demonstrates the real-world application of multilevel models. By estimating the effect of Medicaid expansion while accounting for the nested structure of the data, we can project potential lives saved if non-expansion states adopted the policy. The results show substantial potential benefits, with hundreds of lives potentially saved annually in large non-expansion states.#### 7.2 Identifying High-Risk Areas for Targeted Interventions```{r risk-identification, fig.width=12, fig.height=8}# Identify high-risk counties using model predictions and residualscounty_risk <- county_data %>% mutate( predicted_mortality = fitted(model_cli), state_effect = rep(ranef(model_cli)$state_id[,1], times = counties_per_state), county_residual = residuals(model_cli), risk_score = predicted_mortality + abs(county_residual) ) %>% arrange(desc(risk_score)) %>% mutate(risk_rank = row_number())# Create risk visualizationp_risk_map <- ggplot(county_risk, aes(x = median_income, y = smoking_rate, color = risk_score, size = population)) + geom_point(alpha = 0.6) + scale_color_gradient2(low = "blue", mid = "yellow", high = "red", midpoint = mean(county_risk$risk_score), name = "Risk Score") + scale_size_continuous(name = "Population", range = c(1, 8), labels = scales::comma) + labs(title = "County Risk Profile: Identifying Priority Areas for Intervention", subtitle = "Higher risk scores indicate counties with high mortality and high variability", x = "Median Household Income ($1000s)", y = "Smoking Rate (%)") + theme_minimal() + theme(legend.position = "right")p_risk_map```**Figure Interpretation:** This risk map combines multiple dimensions to identify priority counties for public health interventions. The color gradient (blue to red) indicates overall mortality risk, with red counties having the highest predicted mortality plus unexplained variation. The size of points represents population, helping identify where interventions could have the greatest impact. Counties in the upper-left quadrant (low income, high smoking, shown in warmer colors) represent prime targets for intervention, especially those with large populations. This visualization helps policymakers allocate limited resources to areas with both high need and high potential impact.### 8. Interpreting Multilevel Results and Best Practices```{r interpretation-table}# Create comprehensive interpretation tableinterpretation_data <- data.frame( `Effect Type` = c( "Fixed Effects (County-level)", "Fixed Effects (State-level)", "Cross-level Interactions", "Random Effects (Intercepts)", "Random Effects (Slopes)", "Contextual Effects" ), `What It Measures` = c( "Average relationship between county characteristics and mortality across all states", "Direct effects of state policies and characteristics on county mortality", "How state characteristics modify county-level relationships", "Unobserved state characteristics that create systematic differences in baseline mortality", "How the strength of county-level relationships varies across states", "Difference between within-state and between-state effects" ), `Policy Implications` = c( "Effects of county-level interventions (healthcare access, economic development)", "Effects of state-level policies (Medicaid expansion, health spending, tobacco taxes)", "Whether state policies enhance or diminish county-level interventions", "Need for state-specific baseline adjustments in policy evaluation", "Whether policies should be tailored differently across states", "Whether state context matters beyond individual county characteristics" ), `Example from Our Analysis` = c( "Each $1,000 increase in median income associated with 1.2 fewer deaths per 100,000", "Medicaid expansion associated with 25 fewer deaths per 100,000 on average", "Income effects stronger in non-expansion states (interaction effect)", "States vary in baseline mortality by ±30 deaths per 100,000 due to unmeasured factors", "Income effects vary from -0.8 to -1.6 across states", "State context adds additional protective effect beyond county income" ))interpretation_data %>% kable(caption = "Interpreting Multilevel Model Components in Health Research") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% column_spec(1, bold = TRUE) %>% column_spec(2, width = "25%") %>% column_spec(3, width = "25%") %>% column_spec(4, width = "30%")```### 9. Conclusion and Best PracticesMultilevel modeling provides a powerful framework for analyzing nested health data, offering several key advantages over traditional approaches:#### 9.1. Key Takeaways1. **Hierarchical Structure Recognition**: Health outcomes are naturally nested within geographic, administrative, and organizational units that create dependencies in the data.2. **Variance Partitioning**: Multilevel models partition total variation into meaningful components, helping identify whether interventions should target individuals, communities, or policy levels.3. **Cross-Level Interactions**: These models can test whether higher-level factors (policies, resources) modify the effectiveness of lower-level interventions.4. **Centering Decisions Matter**: The choice between grand-mean centering, group-mean centering, or contextual models fundamentally changes what questions your analysis answers.5. **Policy Evaluation**: The framework supports sophisticated policy evaluation by modeling direct effects, indirect effects, and contextual modifications.#### 9.2. Best Practices for Research```{r best-practices-checklist}best_practices <- data.frame( `Analysis Stage` = c( "Design Phase", "Design Phase", "Design Phase", "Modeling Phase", "Modeling Phase", "Modeling Phase", "Modeling Phase", "Interpretation Phase", "Interpretation Phase", "Interpretation Phase" ), `Best Practice` = c( "Calculate required sample sizes for both levels", "Ensure adequate variation at each level", "Consider balance between levels (e.g., counties per state)", "Start with simple models and build complexity gradually", "Check assumptions thoroughly with diagnostic plots", "Compare multiple model specifications using information criteria", "Consider different centering approaches based on research question", "Focus on policy-relevant effect sizes, not just significance", "Use visualization to communicate multilevel relationships", "Consider substantive significance alongside statistical significance" ), `Health Research Application` = c( "Ensure enough states/regions (>20) and counties/units for reliable estimates", "Verify sufficient variation in policies, outcomes, and predictors", "Balance statistical power with practical constraints", "Build from random intercept → random slope → cross-level interactions", "Examine residuals, random effects normality, and homoscedasticity", "Use AIC, BIC, and likelihood ratio tests for model selection", "Use contextual models when state context theoretically matters", "Report confidence intervals; discuss practical importance for public health", "Create plots showing state-specific relationships and policy effects", "Consider whether effects are large enough to matter for health outcomes" ))best_practices %>% kable(caption = "Best Practices Checklist for Multilevel Health Research") %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% column_spec(1, bold = TRUE, width = "15%") %>% column_spec(2, bold = TRUE, width = "35%") %>% column_spec(3, width = "50%") %>% pack_rows("Study Design", 1, 3, label_row_css = "background-color: #E8F4F8;") %>% pack_rows("Statistical Analysis", 4, 7, label_row_css = "background-color: #F8F4E8;") %>% pack_rows("Results and Communication", 8, 10, label_row_css = "background-color: #F4E8F8;")```#### 9.3. When to Use Multilevel Models in ResearchMultilevel modeling is particularly valuable when:- Data has clear hierarchical structure (counties in states, patients in hospitals)- Research questions involve both individual and contextual effects- Interest lies in understanding variation at multiple levels- Cross-level interactions are theoretically important- Standard errors need to account for clustering- Policy evaluation requires understanding of contextual factors#### 9.4. Common Pitfalls to Avoid1. **Ignoring clustering**: Using standard regression when data is nested leads to incorrect standard errors2. **Over-parameterization**: Adding random slopes for every predictor without theoretical justification3. **Misinterpreting centering**: Not being clear about what centering approach was used and what it means4. **Small sample sizes at Level 2**: Having too few groups (< 20-30) for reliable variance estimates5. **Ignoring diagnostics**: Not checking model assumptions can lead to invalid conclusionsThe combination of multilevel modeling with proper centering strategies, diagnostic checking, and thoughtful interpretation provides health researchers and policymakers with the statistical tools necessary to understand complex, nested relationships in health data and to design more effective, targeted interventions that account for the multilevel nature of health determinants.